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Abstract 


A numerical study of thermohydrodynamic lubrication of a slider bearing has been made. 

The partial differential equations governing the conservation of mass, momentum and en- 
ergy equations coupled with temperature dependent density and viscosity are solved using 
Streamline Upwind Petrov-Galerkin Finite Element Method (SUPGFEM) to yield various 
bearing characteristics. The influence of the solid pad under various thermal boundary con- 
ditions in conjunction with thermal boosting on the load carrying capacity and frictional 
drag force of a slider bearing, at different inlet pressures and for various values of inclination 
parameter have been analyzed. The results depicted through plots and contours indicate 
that a boundary layer and re-circulation zones are noticed in the temperature field. Subse- 
quently the influence of heat conduction both in the stationary pad as well as the moving 
slider on bearing performance have been analyzed. SUPGFEM is effective not only in pre- 
cluding the numerical oscillations arising in the flow field, but also in achieving convergent 
solution. Although various temperature boundary conditions are prescribed on the pad and 
slider, the study reveals that the computed load carrying capacity remains unaffected. This 
is due to the fact that the thermal flux across fluid-solid interface changes marginally only. 
However, when compared with the case of no conduction the load decreases. 

Further, apart from heat conduction in pad, the effect of fluid inertia on bearing char- 
acteristics has also been analyzed. The results reveal that fluid inertia has a profound effect 
on bearing characteristics and its inclusion enhances the load carrying capacity of a slider 
bearing with a reduction in the frictional drag force. 
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Thesis Title: A Numerical Study of Thermohydrodynamic Lubrication of Tilted 
Pad Slider Bearings 

Thesis supervisors: Professor Prawal Sinha and Professor B. V. Rathish Kumar 
Month and year of thesis submission: June, 2000 


The thesis is devoted to the study of thermal and inertia effects on the lubrication char- 
acteristics of a slider bearing. The coupled partial differential equations governing the mass, 
momentum and energy equations of a Newtonian fluid together with temperature dependent 
viscosity and density are solved using the Streamline Upwind Petrov- Galerkin Finite Ele- 
ment Method (SUPGFEM) to yield various bearing characteristics. The influence of thermal 
boosting, various thermal flux boundary conditions, viscosity temperature coefiicient and for 
different values of k (parameter governing inclination of the pad) on load carrying capacity, 
drag force, temperature distribution and on flow fleld have been analyzed. 

The thesis consists of seven Chapters. Chapter 1 gives a General Introduction, wherein 
a brief history of lubrication and various bearing geometries has been presented. It also 
provides a chronological development made in the area of thermal and inertia effects in 
lubrication. This Chapter provides a motivation for the work reported in the subsequent 
Chapters of this thesis. 
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In Chapter 2, the Streamline Upwind Petrov-Galerkin Finite Element Method has been 
discussed. It provides a clear description of Finite Element Numerical Technique and it- 
s applicability for fluid flows. It also presents advantages of FEM in comparison with 
other numerical techniques like Finite Difference Method, Spectral Element Method etc. 
This Chapter also deals with the various weight functions used in the Variational formu- 
lation/Weak formulation and shape functions for the bilinear rectangular elements used in 
SUPGFEM. It also describes the implementation of the SUPGFEM. 

Chapter 3 deals with the investigation of the influence of various temperature boundary 
prescriptions and thermal boosting on bearing characteristics. For the numerical simula- 
tion of hydrodynamic lubrication, we have used Bubnov-Galerkin FEM for a Generalized 
Reynolds equation and SUPGFEM for energy equation in the lubricant film. The study in 
this Chapter reveals that the consideration of (a) Lubricants with low values of viscosity 
temperature coefficient (^5+) (b) Non zero inlet pressure and (c) (inlet temperature) < 

(slider temperature) < (pad temperature) or with non zero flux on the 

pad constitutes a desirable set-up for enhancing the load carrying capacity of a slider bearing 
with a reduction in the frictional drag force. An attempt has also been made for analyzing 
the load generation in a parallel slider bearing (k = 1). 

Lubrication under high load and high speed is common in engineering practice today. 
In the thin film that separates the two rubbing surfaces a large amount of heat may be 
generated. The heat generated in the fluid may escape through conduction in the bearing 
surfaces. In Chapter 4, we have investigated the hydrodynamic lubrication of finite width 
slider bearing taking not only the temperature variation in the fluid film but also the heat 
conduction in the stationary pad. Influence of thermal boosting, different values of inclination 
of the pad and several practical thermal boundary prescriptions on load carrying capacity 
and frictional drag force have been analyzed. SUPGFEM have been used to preclude the 
numerical oscillations arising in the fluid flow. 
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In Chapter 4, we allowed heat conduction only in the stationary surface. But in practice 
heat may escape through the slider also which is moving in the direction of sliding. Hence, 
In Chapter 5, we have included the heat conduction in the pad as well as in the moving 
slider. The system of governing equations which now include two heat conduction equations 
apart from the Generalized Reynolds equation and energy equation are solved simultaneously 
using SUPGFEM. The results indicate that this method is effective not only in precluding the 
numerical oscillations in the flow field but also obtaining a convergent solution. Load carrying 
capacity of slider bearing can be thermally boosted by setting slider at temperature lower 
than that of incoming lubricant. In this study we have chosen different boundary settings 
for pad and slider conduction. These boundary conditions do alter the temperature field in 
the pad. However, the thermal flux across the fluid-solid interface changes marginally only. 
Consequently the load carrying capacity of the slider bearing remains unaffected. It is to be 
noted that the isothermal pad employed in the present study has been set at temperature 
similar to that of inlet flow. 

In Chapter 6, we have included the effect of fluid inertia along with the thermal effects. 
The system of nonlinear partial differential equations governing the conservation of mass, 
momentum and energy with viscosity and density depending on the temperature has been 
solved. Prom the results, it is observed that fluid inertial effects are significant when the 
Reynolds number exceeds 5 x 10^. Inclusion of fluid inertia enhances the load capacity of a 
slider bearing for various values of k and at various inlet pressures. Also, the investigation 
on the influence of viscosity temperature coeflflcient {P'^) on the load carrying capacity at 
various inlet pressures has been made. The dependency of load carrying capacity on for 
pf = 0, 0.5, 1 shows that, for small inlet pressures, load generation in slider bearing decreases 
with increasing but at high inlet pressures this reduction in load generation becomes 
marginal. It can also be noted that for fixed T'*', increases with decreasing p'^ and 

thus the effective viscosity of the lubricant increases and this leads to a higher load capacity. 
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Chapter 7 deals with the combined influence of fluid inertia and heat conduction in the 
pad. In this Chapter, load carrying capacity for various thermal boundary conditions are 
shown in the form of Tables. The computed load carrying capacity have been compared 
with the results obtained for no conduction in the pad. The results show that for higher tilt 
(decreasing k) load capacities are high for the case with no conduction. Whereas for lower 
tilt (increasing k) the load carrying capacities are low for the case with conduction. This 
could be due to an increase in the lubricant temperature which may not have been able to 
escape as the conduction in the pad was not allowed. Also the nonlinear model (with inertia) 
has an upper hand as compared with the linearized model (with no inertia) in the computed 
load capacities for various values of k. 
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Chapter 1 

General Introduction 


1.1 Lubrication 

Lubrication is the science and technology of interacting surfaces in relative motion and of 
the practices thereto. It is simply the use of a material to improve the smoothness of 
movement of one surface over another surface. The material which is used to achieve this 
is called lubricant and the process is called “Lubrication.” The most common lubricants 
are oils and greases but other kinds of liquid substances as well as solids and gases are also 
used as lubricants. Among the solid lubricants graphite, borax, mica, white lead, waxes 
etc. are found to be used in deep drawing, extrusion and light rolling operations in metal 
industries, and many modern engineering mechanics, especially in the field of nuclear power, 
aviation and ballistic missiles. Gases are finding increasing usage in gas-cycle machinery, 
gyros, where precision and constancy of torque are critical; in food and textile processing 
machinery, where cleanliness and absence of contaminants are critical; and in high speed 
dental drills. Lubricants in general help in reducing or controlling friction. However a 
reduction is not always desirable. There may be situations in which it is more important 
to maintain steady friction than to obtain the lowest possible friction. Some examples are, 
the control of chatter in a machine tool slide ways or grinding operation, control of strip 
movement in metal rolling and elimination of brake sequel or clutch judder in a car. In 
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addition to controlling friction, lubricants are usually expected to reduce wear and often to 
prevent overheating and corrosion. In some cases, lubricants are also used as a cooling agent 
by carrying away heat which is generated in bearings. For example, the chief requirement of 
the lubricants in forging are cooling the dies and lubrication of the die-workpiece interface 
in order to minimize metal pick-up, reduce wear and facilitate plastic flow of metal. 


1.1.1 Historical Perspective of Lubrication 


The understanding of the hydrodynamic lubrication began with the classical experiment of 
Tower [1] , in which the existence of a film was detected from measurements of pressure 
within the lubricant, and of Petroff [2, 3, 4, 5] , who reached the same conclusion from 
frictional measurements. Tower [6] demonstrated that a loaded, oil lubricated journal bearing 
is subject to local pressures substantially higher than its mean pressure. Tower’s [1, 6] 
experiments to study friction in railroad journal bearings, eventually led to the development 
of the theory of hydrodynamic lubrication. The work was closely followed by Reynolds [7] and 
he used linearized form of the equations of motion for laminar flow of a Newtonian fluid with 
constant viscosity, between rigid surfaces and derived a second order differential equation for 
the pressure in the narrow converging gap between bearing surfaces. This pressure enables 
a load to be transmitted between the surfaces with extremely low friction, since the surfaces 
are completely separated by a fluid film. In such a situation the physical properties of the 
lubricant, notably the dynamic viscosity, dictate the behavior. The equation describing the 


pressure in two dimensional hydrodynamic bearings is given by 


d .h^ dp. 
dx^ 7] dx 


6Uo 


dx 


( 1 . 1 ) 


where x is the coordinate in the sliding direction, p is the pressure, rj is the viscosity, h is the 
film thickness and C/q is the effective sliding velocity. The work of Reynolds received imme- 
diate acceptance by the engineering community, but there were great difficulties encountered 
when attempting to solve the Reynolds equation. As numerical solutions of the Reynolds 
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equation became a common place, more and more investigators began to ask the question 
i.e. Does the Reynolds pressure equation account for the physics of the problem adequately?. 
In fact the extension of classical theory was found to be wanting much earlier. Martin [8] 
found that it cannot satisfactorily explain the existence of continuous hydrodynamic film 
in highly stressed contacts, such as those found in gears. With the increasing demands of 
various sophisticated engineering equipments containing bearings, such as metal industries, 
die industries and many other application, the extension of Reynolds classical theory became 
necessary, to predict the bearing operations accurately. 

1.1.2 Bearings and their Characteristics 

A bearing is a support or a guide that locates one machine component with respect to others 
in such a way that prescribed relative motion can occur while the forces associated with 
machine operation are transmitted smoothly and efficiently. Bearings can be classified in 
several ways: according to the basic mode of operation (rubbing, hydrodynamic, hydrostatic, 
or rolling element), according to the direction and nature of the applied load (thrust or 
journal), or according to geometric form (tapered or stepped parallel surface or tilted pad). 
Various types of bearings are shown in Figures 1.1 - 1.5. 


o ol= e (eccentricity) 

N = shaft speed 
(R-r) = Cr 
n = e / Cr 

h = Cr ( 1 + n Cos 9 ) 



Figure 1.1: Journal bearing 
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Figure 1.2: Parallel surface slider bearing 




Figure 1.4; f'ixed-incline (tilted pad) slider bearing 
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P 



Figure 1.5; Parallel step slider bearing 

The basic function of any bearing, irrespective of its configuration or working principle 
is to restrain unwanted relative movement between the two machine parts while offering the 
least resistance to the desired relative movement. The other characteristics of a bearing 
include: 

• Load carrying capacity 

• Stiffness 

• Flow rate of supporting fluid 

• Resistance to sliding 

• Power required to operate 

• Heating. 

Design is a creative process aimed at finding a solution to a particular problem. In all forms 
of design a particular problem may have, many different solutions, mainly because design 
requirements can be interpreted in many ways. For example, it may be desirable to produce: 

• The cheapest solution 

• The bearing with ease and with variable material 
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• The most reliable 

• The one that takes up the smallest space 

• The one that is lightest in weight 

• The best from any of a whole variety of possible stand points. 

Therefore, the process of selection and design usually involves these steps. 

• Selecting a suitable type of bearing 

• Estimating a bearing size that is likely to be satisfactory 

• Analyzing bearing performance to see if it meets the requirements 

• And then modifying the design and the dimensions until the performance is near to 
whichever optimum is considered the most important. 

1.1.3 Lubrication Regimes 

In general the bearings operate under one of the following modes of lubrication. 

• Hydrodynamic/Hydrostatic lubrication 

• Boundary lubrication '* 

• Mixed lubrication. 

In hydrodynamic lubrication, moving surface(s) of the bearing are completely separated by 
a layer of fluid film of thickness of the order 10“’'m - 10“®m and the pressure in the film is 
generated by hydrodynamic action. Such bearings are called hydrodynamically lubricated 
bearings. In hydrodynamic- lubrication the films are generally thick so that opposing solid 
surfaces are prevented from coming into contact. This condition is often referred to as “The 
ideal form of lubrication”, since it provides low friction and high resistance to wear. The 
lubrication of the solid surfaces is governed by the bulk physical properties of the lubricant, 
notably the viscosity, and the frictional characteristics arise purely from the shearing of the 
viscous lubricant. For a normal load to be supported by a bearing, positive-pressure profiles 
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must be developed over the bearing length. For a positive pressure to be developed in a slider 
bearing (see Figure 1.6 (a)) the lubricant film thickness must be decreasing in the sliding 
direction. In a squeeze film bearing (see Figure 1.6 (b)) the bearing surfaces approach each 
other with a normal velocity (squeeze velocity) Ua- The squeeze mechanism of pressure 
generation provides a valuable cushioning effect when the bearing surfaces approach each 
other. Positive pressures will be generated only when the film thickness is diminishing. 

In some situations, load on the bearing is too high and the relative velocity of the 
bearing surfaces is too low to generate sufficient hydrodynamic pressure in the oil film. In 
such cases, an external pressure source may be used to supply oil under a high pressure to 
prevent metal - to - metal contact between bearing surfaces. Such mode of lubrication is 
called externally pressurized or hydrostatic. In an externally pressurized bearing (see Figure 
1.6 (c)) the pressure across the bearing supports the load. 



(C) 

Figure 1.6; Mechanism of pressure development for hydrodynamic lubrication (a) Slider bear- 
ing (b) Squeeze film bearing (c) Externally pressurized bearing 


In case of boundary lubrication the solid surfaces are not separated by a thick lubricant 
film, fluid film effects are negligible and there is considerable asperity contact. The frictional 
characteristics are determined by the properties of the solids and the lubricant film at the 
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common interfaces. The surface films vary in thickness from Inm to lOnm depending 
on the molecular size. In boundary lubrication, although the friction is much higher than 
in the hydrodynamic regime, it is still much lower than for unlubricated surfaces. The 
contact mechanism is governed by the physical and chemical properties of thin surface films 
of molecular proportions. The properties of the bulk lubricant are of minor importance, and 
the friction coefficient is essentially independent of fluid viscosity. Figure 1.7 illustrates the 
film conditions in fluid film and boundary lubrication. The surface asperities in this Figure 
are greatly distorted for purposes of illustration. To scale, real surfaces would appear as 
gently rolling hills rather than the sharp peaks. The surface asperities are not in contact for 
fluid film lubrication but are in contact for boundary lubrication. 

mrfYrm 

(a) 

Figure 1.7: Film conditions of lubrication regime (a) Fluid film lubricated-surfaces separated 
by bulk lubrication film (b) Partial lubrication-both bulk and boundary lubrication play a role 
(c) Boundary lubrication-performance depends essentially on boundary film 

Mixed lubrication is a situation between hydrodynamic and boundary lubrication. The 
fluid film between the two moving surfaces is thin enough and the surface asperities begin 
to interfere in the hydrodynamic process. The load is partly supported by asperity contacts 
and partly by fluid film. Figures 1.8 and 1.9 shows the friction coefficient and wear rate 
respectively for various types of lubrication regimes. 
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Figure 1.8; Bar diagram showing friction coefficient for various lubrication regimes 



Figure 1.9: Wear rate for various lubricating regimes 
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1.2 Extension of Classical Lubrication Theory 

There were at least six areas in which the extension of classical theory seemed desirable. Each 
of these areas is concerned with the relaxation of one or more constraints of the classical 
theory and leads to studies of: 

• Thermohydrodynamic effects 

• Inertia effects 

• Compressibility effects 

• Elastohydrodynamic effects 

• Turbulence effects 

• Non-Newtonian effects. 

This thesis is mainly concerned with the first two effects. Hence, in the subsequent Chapters 
emphasis is given to the related investigations only. 

1.2.1 Thermohydrodynamic Effects 

In recent years, there has been a growing demand for bearings which can sustain higher load 
and are at relatively higher speeds. A large amount of heat is generated in such bearings. 
This necessitates the need for fundamental understanding of the thermal conditions in the 
bearings. Under normal operating conditions most lubricants undergo an increase in tem- 
perature due to viscous dissipation which in turn may significantly influence the viscosity 
and load carrying capacity of the bearing. Two parameters which are generally responsible 
for bearing failure are 

• Maximum bearing temperature 

• Minimum film thickness. 
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In view of this, the actual temperature distribution at every point needs to be estimated. This 
can be obtained through the simultaneous solution of equations of continuity, momentum and 
energy conservation. One of the difficulties encountered in dealing with energy conservation 
stems from the fact that there are no exact analjdical expressions for the viscosity as a 
function of temperature and pressure. This is because viscosity variation occurs mainly due 
to 


• Heat generation by viscous shear 

• Heat transfer to and from the surface. 

In many applications, increase in temperature across the bearing, although not negligible, 
remains small. The classical isothermal theory in this case might still be applied as a first 
approximation, provided that calculations are based on an “effective” or “operating” vis- 
cosity. In large bearings and/or under severe running conditions, thermal effects may be 
si gnifi cant (ref: Pinkus, [31], page; 4) and prediction of bearing performance may no longer 
be possible when based on the assumption of a uniform effective viscosity. Thus, numerical 
emperical relations of /x = /x(T) has been developed and the one, which is extensively used 
is given by /x = /xo exp [-^(T - To)]; where is the viscosity at temperature To and P is 
a constant, called the viscosity temperature exponent. In view of these aspects, the thermo- 
hydrodynamic model can be defined as “A model which predicts bearing performance based 
on realistic thermal boundary conditions on the stationary and moving surfaces and includes 
a pointwise variation of the fluid viscosity with temperature” . The study of thermohydro- 
dynamic model involves simultaneous solution of energy balance, the conservation equations 
of mass and momentum which can be achieved by an iterative numerical process. 

The early applications of energy equation to lubrication theory were by Christopher- 
son [9] and by Cameron and Wood [10] . The hydrodynamic theory of film lubrication was 
investigated by Cope [11], who systematically derived the governing equations of fluid film 
lubrication with thermal effects. Assuming that oil carries away the entire frictional heat, he 
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derived a simplified energy equation in which the variation of viscosity across the film was 
neglected. Hunter and Zienkiewicz [12] showed that in certain circumstances, the viscosity 
variation across the lubricant film was responsible for generation of an applicable load and, 
therefore, should not be neglected. Dowson [13] presented a Generalized Reynolds equa- 
tion and considered the energy equation by taking the variation of viscosity both along and 
across the lubricant film. Several other researchers who studied the thermal effects include 
Huebner [14], Ezzat and Rhode [15], Khonsari [16], Gethin [17], Gero and Ettles [18], 
Schumack [19] etc. 

A review of major contributions on thermal effects in various geometries of bearings is done 
by Khonsari [16]. He has presented an extensive literature survey on heat effects in various 
type of bearings. Summary, given by Khonsari regarding the advances made in analyzing 
the thermal effects in bearings, conclude: 

• Thermal effects are important as they influence bearing perforihance significantly 

• Viscosity variation as a function of temperature must be considered in bearing design 

• The effect of temperature on viscosity is more pronounced than the effect of elastic 
distortions due to temperature and pressure 

• Under most operating conditions, the heat conduction along the direction of motion is 
small and can be neglected. However, when back flow of oil, is produced the conduction 
term along the direction of motion must also be retained in the energy equation 

• Cooling the slider is generally advantageous because it can increase the load carrying 
capacity and reduces the friction which is one of the most interesting aspect for the 
designers. 

1.2.2 Inertia Effects 

The performance of a typical bearing is significantly affected by inertia effects that are 
generated w'ithin the lubricant. This complex mechanism stems from the convective nature 
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of the flow fleld, and are driven by a combination of high operating speeds, low viscosity 
of the lubricant, large external pressures, and the complexity of the bearing geometry. For 
the bearing recess region, convection effects will lead to a non-linear pressure distribution. 
It was found that in typical hydrodynamically lubricated thrust and journal bearings the 
inertia effects are negligible as compared to the viscous effects. However, inertia effects 
can sometimes be important. Two examples are given in Myllerup and Hamrock [20], the 
first example is for any surface feature resulting in a film geometry where || Order 
(1) or ^ PS Order (1) where h = h{x) is film thickness. The singularity arising from 
the curvature is less severe than the one arising from the film gradient, but in both cases 
using the Reynolds equation will not allow a proper matching between the fluid regions. The 
second example where the Reynolds equation will not describe what occurs in the lubricant 
conjunction is when a solid particle is within the conjunction. In these cases, it is necessary 
to consider the convective inertia. Inertia effects were studied by Tichy [21] , Talmage et 
al. [22], Prasad et al. [23] etc. All these studies concluded that inertia effects are more 
significant when the Reynolds number is not small and bearing performance is influenced by 
the mass effect of the lubricant and high operating speeds. 

1.2.3 Compressibility Effects 

When the fluid (lubricant) is compressible (e.g. Gases), compressibility effects have to be 
considered. Gas lubricated film bearings have commanded considerable attention in recent 
years because they possess characters that make their usage advantageous in many bearing 
application. They are analogous to the hydrodynamic oil-lubricated bearings except that the 
fluid is compressible. Further more, since air is 1000 times less viscous than even the thinnest 
mineral oils, the viscous resistance is very much less. However, the distance of nearest 
approach between the bearing surfaces is correspondingly smaller so that special precautions 
must be taken in manufacturing the bearing. Gross [24], Houpert and Hamrock [25] etc. 
have given an extensive study about the compressibility effects. 
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1.2.4 Elastohydrodynamic Effects 

In recent years a number of practical bearings have employed layered elastic solids as mating 
surfaces. Also to avoid wear and seizure and to make the supply of the lubricant easier, thin 
layer of softer material, such as stilfer plastics, are sometimes coated on the bearing surface. 
In such situations, neglecting the surface deformation will not predict the film profile and 
pressure distribution accurately. This has led to the substantial interest in the elastohy- 
drodynamic lubrication of such bearings and it is imperative to consider the deformation of 
bearing solids. 

When the surfaces are no longer rigid, elastohydrodynamic effects become significant. 
The viscosity must be considered to be pressure dependent. Elastohydrodynamic lubrication 
can be defined as the study of lubricating film between elastically deformable solids. At one 
end of the scale, it covers highly stressed lubricated machine elements like ball and roller 
bearings and gears and at the other it includes the performance of soft rubber seals and 
bearings with soft liners. Hamrock and Dowson [26, 27, 28] Ghosh and Hamrock [29] etc. 
have given extensive explanation on elastohydrodynamic lubrication. 

1.2.5 Turbulence Effects 

It is well known that the motion of fluids may occur in ‘irregular fluctuations’, resulting 
in whirls and/or eddying. Such motions are called turbulent flows. When the Reynolds 
number {Rg = is increased (outside its range for laminar flow), the boundary layers 
around solid bodies change from laminar to turbulent. Such a transition is influenced by 
geometries, pressure gradients, suction, compressibility, and heat transfer. It is mandatory 
to consider fluid inertia along with the turbulence. This would indeed carry the governing 
equations for fluid film lubrication beyond the customary postulates of lubrication theory. 
Extensive study have been considered by Venkateswarlu and Rodkiewicz [30], Pinkus [31] 


etc.- 
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1.2.6 Non-Newtonian Effects 

For most fluid film lubrication analysis the lubricant is assumed to behave in a Newtonian 
manner. This implies that the shear stress r is linearly related to the shear strain s. i.e 
s =z I. Where rj is the absolute viscosity or proportionality constant. It may vary with 
both pressure and temperature. However, in many situations the lubricant may experience 
rapid and extremely large pressure variations, large temperature changes, and particularly 
in sliding contacts, high shear rates, this implies that the lubrication shear stress is still a 
function of the shear rate but that the relationship may no longer be linear. Such class of 
fluids are called non-newtonian fluids. Shukla and Isa [32] , Prasad et al. [23] have given 
an extensive study. Hirst and Moore [33] presented an experimental investigation of the 
dependence of shear stress and shear strain rate and proposed a Eyring type of relationship. 
Several other non-notonian models have been proposed, such as Power law fluids (Prasad et 
al. [23]), Bingham fluids (Caldewell and Babbitt [34]), Viscoplastic fluids (Sheu et al. [35]) 
etc. 


1.3 Review of the Literature 

It was Fogg [36] who first published his experimental findings suggesting that a pressure film 
is generated, in a slider bearing, as a result of thermal expansion of the fluid passage through 
the bearing. Subsequently, this phenomena received considerable attention under various 
titles, like “thermal wedge” , “viscosity wedge” , or “density wedge” . In the Fogg experiment 
radial grooves were made as sharp as possible in the direction of rotation of the moving collar 
and the direction of rotation of moving collar could be reversed. It was concluded that a 
parallel surface without a rounded inlet edge was capable of sustaining considerable loads 
which was attributed to the thermal expansion of the fluid. This appears to be plausible 
except that Fogg [36] recognized only one contributing factor while there must have been, 
at least, two factors generating lifting force. The other factor is the pressure build-up in 
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the fore-region (Rodkiewicz et al. [37]). This must have resulted in an above ambient inlet 
pressure, leading to a non-zero load carrying capacity. 

Osterle et al. [38] specialized their governing equations to the infinitely wide parallel- 
surface slider bearing (k = 1), according to presentation of Cope [11]. Inertia forces were 
neglected, pressure was assumed constant across the fluid film thickness and viscous force 
was represented by the dominant term in the momentum equation. However, the density 
and viscosity were not assumed to vary exponentially with temperature and in the energy 
equation the conduction terms were disregarded. Lewicki [39] introduced the concept of 
pressure build-up in the inflow zone of a square-faced slider. In his incompressible flow 
analysis with constant viscosity the stream function approach is used within the inflow zone, 
and within the lubricating film, pressure is assumed to vary linearly. 

The effect of prescribed thermal gradients was also studied for journal bearings as well as 
for slider bearings by Hagg [40], Tipei and Nica [41]. Zienkiewicz [42] analyzed an infinitely 
long parallel slider bearing allowing for temperature variations across the film. Cameron and 
Wood [10] considered “the thermal wedge” , and Cameron [43] “viscosity wedge” in a thrust 
bearing. Basically their problem formulation is similar to that of Osterle et al. [38]. They 
also disregard the conduction terms in the energy equation. Hunter and Zienkiewicz [12] 
performed analysis for the one-dimensional slider bearings. In this analysis the thermal 
boundary conditions at the fluid-solid interface were assumed a priori. Dowson and Hud- 
son [44] re-examined the same problem allowing for the heat transfer to the bearing solids. 
The effect of thermal distortions of the bearing solids was considered by Hahn and Ket- 
tleborough [45] in an analysis of the one dimensional slider bearing. They concluded that 
while thermal distortions are responsible for the load carrying capacity exhibited by parallel 
surface slider; they have little effect on the performance of inclined sliders. Also Hahn and 
Kettleborough [46] presented an iterative method of calculating the steady-state pressure and 
temperature distributions for infinitely wide slider bearings. The viscosity was a function of 
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temperature and pressure, and the density a function of temperature. 

Lebeck [47, 48] have experimentally investigated the thermal influence on load carrying 
capacity of parallel slider bearings. Rohde and Kong [49] investigated the thermal influence 
on bearing characteristics using the finite difference numerical technique. The effect of fluid 
inertia on bearing characteristics has been analyzed by Hill et al. [50], Artiles et al. [51]. In 
these investigations, it is concluded that inclusion of inertia has a profound effect on bearing 
characteristics. Yang and Rodkiewicz [52] also investigated the influence of thermal effects 
in tilted pad slider bearings. 

Another important factor which influences bearing characteristics in case of high speed 
bearing lubricated with low viscosity fluid and subjected to a squeezing motion is the fluid 
inertia (Hashimoto and Wada [53], Elkouh [54]). This effect is important even in case of 
laminar flow where an increase in hydrodynamic film thickness may also some times compel 
one to consider this effect. Dowson et al. [55] conducted a wide range of experiments on a 
cylindrical-plane bearing. It was observed that for a moderately large Reynolds number, the 
experimental data deviated substantially (due to fluid inertia) from the theoretical values 
predicted by classical Reynolds equation using the Coyne-Elrod film rupture conditions. 
You and Lu [56] presented a theoretical analysis of a cylindrical-plane bearing and journal 
bearing. They have shown that the lubricant inertia has a profound effect even at moderate 
values of the Reynolds number. 

Some recent work in this direction is by Talmage et al. [22], Zhang and Rodkiewicz 
[57] who considered the influence of fluid inertia on bearing characteristics using a CFD 
technique called Pseudo spectral finite difference method. Their work was mainly on scheme 
development which could simultaneously solve the complete Navier-Stokes equations, the 
continuity equation and the energy equation. This method can accommodate arbitrary inlet 
and exit boundary conditions as in a coupled analysis. But the approach requires greater 
effort and computational resources. And in all these studies the compressibility work term 
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was ignored. 

Recently, Rodkiewicz and Sinha [58], Schumack [19] have numerically analyzed the 
problem of thermohydrodynamic lubrication with Dirichlet boundary prescription on the pad 
and slider. But the actual slider bearing configuration would require zero or non zero flux 
boundary prescriptions. Gero and Ettles [18] has attempted the same problem but his study 
was focussed on scheme development. Further, he has dealt with either constant viscosity 
or density or linearized viscosity with constant density. The influence of thermal boosting, 
which amounts to considering slider, pad and inflow lubricant at different temperatures, on 
load carrying capacity, drag force etc. of the lubricant in thermohydro dynamic slider bearing 
context has not been investigated. 

In view of the above comments, a thorough investigation has been made in this thesis to 
analyze the bearing characteristics in a slider bearing lubrication. Also particular emphasis 
has been made for the parallel surface slider bearing lubrication. 

1.4 Outline of the Dissertation 

The main objective of the present study is to simulate the load carrying mechanism in slider 
bearings using the Streamline Upwind Petrov-Galerkin Finite Element Method (SUPGFEM), 
because the applicability of FEM alone is limited in cases where the governing equations are 
highly non-linear, such as those found in the lubrication theory with thermal effects. 

The thesis consists of seven Chapters, where Chapter 1 is the General Introduction, 
which provides a brief history of lubrication and outline of various lubricating regimes. It 
also presents a chronological developments made in the area of thermal and inertia effects. 
This Chapter highlights the motivation for the work reported in the subsequent Chapters. 

In Chapter 2, we have discussed various forms of obtaining a approximate numerical 
solution for a differential equation. Special attention has been made for the methodology 
“Streamline Upwind Petrov-Galerkin Finite Element Method”, its applicability for fluid 
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flows and various steps involved in the implementation of the method has been discussed 
extensively. 

Chapter 3 describes the mechanism of tilted pad slider bearings and influence of various 
temperature boundary conditions on bearing characteristics. It also deals with the influence 
of various inlet pressures, different values of k, /?+, various slider and pad temperature 
settings on the velocity distribution, temperature distribution, load carrying capacity and 
frictional drag force. Thorough discussion of iso-u velocities and isotherms have also been 
provided. 

In Chapter 4, we have investigated the hydrodynamic lubrication of finite width slider 
bearing considering not only the temperature variation in the fluid film but also the conduc- 
tion. Influence of thermal boosting for various values of k and several practical temperature 
boundary conditions on load carrying capacity and frictional drag force has been discussed. 
SUPGFEM is found to be very effective in precluding the numerical oscillations arising in 
the fluid flow. 

In Chapter 4, we have included only heat conduction in the stationary pad. But in 
practice the heat generated in the fluid may escape either through conduction in the pad 
or slider or both. So in Chapter 5, we have analyzed the performance of slider bearing 
considering heat conduction in both the bearing surfaces (slider and pad). SUPG weight 
functions are found to be advantageous in obtaining a convergent solution. In this Chapter, 
pad and slider are set at various boundary conditions. These boundary conditions do not 
have a significant variation in their respective computed load capacity, because the variation 
in the temperature flux gradients is negligibly small. 

In Chapter 6, the influence of convective inertia has been included. The nonlinear 
partial differential equations governing the conservation of mass, momentum and energy with 
temperature dependent density and viscosity are solved using SUPGFEM to yield various 
bearing characteristics. Also we have investigated the influence of viscosity temperature 
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coefficient on load carrying capacity. In this Chapter, we have discussed iso-u velocity 
contours and isotherms for various values of T^,T^,T4^. 

Chapter 7, deals with the combined influence of convective inertia and heat conduction 
in the stationary pad. In this Chapter, a comparison has been made for the efiPective load 
carrying capacity values with conduction in the pad and no conduction in the pad. 


Chapter 2 

Streamline Upwind Petrov- Galer kin 
Finite Element Method 


2.1 Introduction 

While searching for a quantitative description .of physical phenomena, the engineer or the 
physicist establishes generally a system of ordinary or partial differential equations valid in 
a certain region (or domain) and imposes on this system, suitable boundary and/or initial 
conditions. At this stage the mathematical model is complete, and for practical applications 
“merely” a solution for a particular set of numerical data is needed. Here, however, come 
the major difficulties, as only the very simplest forms of equations, within geometrically 
trial boundaries, are capable of being solved exactly with available mathematical methods. 
Ordinary differential equations with constant coefficients are one of the few examples for 
which standard solution procedures are available and even here, with a large number of 
dependent variables, considerable difficulties are encountered. 

To overcome such difficulties and to enlist the aid of the most powerful tools developed 
on computer based techniques, it is necessary to recast the problem in a purely algebraic 
form, involving only the basic arithmetic operations. To achieve this, various form of the 
continuum problem defined by the differential equations can be used. In such a discretization 
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the infinite set of numbers representing the unknown function or functions is replaced by a 
finite number of unknown parameters, and this process, in general, requires some form of 
approximation. 

Of the various forms of discretization which are possible, one of the simplest is the 
finite difference process. In the finite difference approximation of a differential equation, the 
derivatives are replaced by difference quotients (or the function is expanded in a Taylor’s 
series) that involve the value of the solution at discrete mesh points of the domain. The 
resulting algebraic equations are solved, after imposing the boundary conditions, for the 
values of the solution at the mesh points. So the finite difference model of a problem gives 
a pointwise approximation to the governing equations. The model is improved (formed by 
writing difference equations for an array of grid points) as more points are used. With finite 
difference technique, we can treat some fairly difficult problems; but, for example, when we 
encounter irregular geometries or unequal specification of boundary conditions, we find that 
finite difference technique becomes hard to use. And also it is not economical either in terms 
of computer storage or CPU time. 

For instance, in the solution process of the incompressible lubrication problems for 
bearings of complex geometry, the finite difference method is inherently hampered by many 
limitations. In the technique for each problem of interest, a preferred coordinate system is 
found so that lubricant boundaries are described by constant coordinate lines, or a loss in 
accuracy of the solution is accepted. In addition, when abrupt variations of the field prop- 
erties, such as film thickness, occur, auxiliary conditions, such as continuity, are generally 
invoked. 

Thus, the need for a general yet flexible method that will adequately treat diverse 
problems, within a single framework of numerical procedures is clearly apparent. One such 
procedure is finite element method (FEM). Unlike the finite difference method, which envi- 
sions the solution region as an array of grid points, the finite element method envisions the 
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solution region as built up of many small, interconnected subregions or elements. A finite 
element model of a problem gives a piecewise approximation to the governing equations. 
The basic premises of the finite element method is that a solution region can be analytically 
modeled or approximated by replacing it with an assemblage of discrete elements. Since 
these elements can be put together in a variety of ways, they can be used to represent ex- 
ceedingly complex shapes. As an example of how finite difference and finite element model 
might be used to represent a complex geometric shape. Consider the turbine blade cross 
section shown in Figure 2.1. 




Figure 2.1: Finite difference and finite element discretization of a turbine blade profile (a) Typical 
finite difference model, (b) Typical finite element model. 

For this device, we may want to find the distribution of displacements and stresses for 
a given force loading. A uniform finite difference mesh would reasonably cover the blade 
(the solution region), but the boundaries must be approximated by a series of horizontal 
and vertical lines. On the other hand, the finite element model (using the simplest two- 
dimensional-the triangle) gives a better approximation to the region and requires fewer nodes. 
Also, a better approximation to the boundary shape results because the curved boundary is 
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represented by a series of straight lines. And the grid need not be structured. Due to this 
unstructured form very complex geometries can be handled with ease. This is clearly the 
most important advantage of finite element method and is not shared by the familiar finite 
diflFerence method, which needs a structured grid. 

The first essential characteristic of the finite element method (FEM) is that the solution 
of the discrete problem is assumed a priori to have a prescribed form. The solution has to 
be in a function space, which is built by varying function values in a given way, for instance 
linearly or quadratically, between values in nodal points. The nodal points or nodes, are 
typical points of the elements such as vertices, mid side points or mid-element points, etc. 
Due to this choice, the representation of the solution is strongly linked to the geometric 
representation of the domain. 

The second essential characteristic is that FEM does not look for a solution of the partial 
differential equation (PDE) itself, but looks for a solution of some integral form of PDE. In the 
most general form, this integral formulation is obtained from a weighted residual formulation. 
By this formulation the method acquires the ability to naturally incorporate differential type 
boundary conditions. This property constitutes the most important advantage of the FEM 
and is not shared by any other method (e.g. Finite Difference Method, Finite Volume 
Method, Spectral Element Method) . 

The combination of the representation of the solution in a given function space, with 
the integral formulation treating rigorously the boundary conditions, gives to the method an 
extremely strong and rigorous mathematical foundation. This allows, for instance, a precise 
definition of accuracy. The concept of accuracy is very loosely defined in FDM where it often 
depends upon strongly on the analytic regularity of the solution. 

The third essential characteristic of the FEM is the modular way in which the discretiza- 
tion is obtained. The discrete equations are constructed from the contribution on element 
level which afterwards are assembled to get an over all study of a particular problem. 
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2.2 Finite Element Method 

The finite element method is a type of a scheme which is used to reduce a field problem 
described by an extremum principle to one in which there are a finite number of unknowns. 
The method uses many localized functions, termed interpolation functions, each valid only 
for a small subregion of the domain of interest. Thus the overall state of the unknown field 
variable is represented piecewise, with continuity of the variable across the boundaries of the 
subregions or elements being assured by suitable requirements placed on the choice of the 
interpolation polynomials. Clearly, discretizing the region of the field problem into a smaller 
regions does not by itself reduce the problem as the continuum to a problem of a finite 
number of unknowns; within each subregion a continuum problem still exists. However, this 
is where the finite element method makes the most significant assumption that leads to its 
success in solving continuum problems. It assumes that the state of the field variable within 
the subregion or element is completely described by values of the unknown variable at a 
finite number of points. 

The method is a valuable tool in the solution of many engineering problems. In situations 
where the governing equations are known, but complicated geometry or boundary conditions 
render analytical solutions difficult or impossible to obtain, the finite element method is often 
employed. The finite element method makes use of spatial discretization and a weighted 
residual formulation to arrive at a system of matrix equations. Solution of the matrix 
equations yields an approximate solution to the original boundary value problem. The most 
common weighted residual formulation has been the Galerkin method, in which weighting 
and interpolation functions are from the same class of functions. The Galerkin method, when 
applied to most structures or heat conduction problems, leads to symmetric stiffness matrices. 
In this case, it can be shown that the solution possesses ‘the best approximation’ property. 
That is, the difference between the finite element and the exact solution is minimized with 
respect to a certain norm. The success of the Galerkin finite element method in structural 
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applications is largely due to the ‘best approximation’ result. 

Compared to other numerical techniques, the finite element method possesses numerous 
advantages. In addition to complete generality with regard to geometry and field property 
variations, prescription of boundary conditions in terms of pressure or flow, cyclic conditions, 
present no special difficulties. In the limit, the method can be shown to converge to the 
exact solution; Thus any degree of accuracy can be obtained. Moreover, the system of 
equations generated by the method produces a symmetric, banded, definite matrix which 
can be solved directly with minimal requirements of computer storage and time. The method 
is characterized by three features. 

• The domain of the problem is represented by a collection of simple subdomains, called 
finite elements. The collection of finite elements is called the finite element mesh 

• Over each finite element, the physical process is approximated by functions of desired 
type (polynomials or otherwise), and algebraic equations relating physical quantities 
at selective points, called nodes of the element are developed 

• The element equations are assembled using continuity and/or “balance” of physical 
quantities. 

In the finite element method, in general, we seek an approximate solution u to a differential 
equation in the form 

n 

w ~ E %% 

where Uj are the values of u at the element nodes, iJjj are the interpolation functions. Direct 
substitution of such approximation into the governing differential equations may not always 
result, (for an arbitrary choice of the data of the problem) in a necessary and sufficient 
number of equations for the undetermined coefficients uj. Therefore, a procedure where by a 
necessary and sufficient number of equations can be obtained is needed. One such procedure 
is provided by a weighted-integral form of the governing differential equation. 




2.3 Variational Methods 


27 


2.2.1 Historical Perspective of FEM 

Historically, the finite element method originates from structural mechanics. The first effort 
to use FEM via piecewise continuous function defined over triangular domains appear in 
the applied mathematics literature with the work of Courant [59]. Motivated by Euler [60], 
Courant used an assembly of triangular elements and the principle of minimum potential 
energy to study the St. Vennat torsion problem. After Courant’s work, nearly a decade 
passed before these discretization ideas were used again. The work of Weinberger [61, 62], 
Poly [63, 64], Hersch [65] who focussed his attention on boundary eigen values, were a period 
of renewed interest. 

As the popularity of the FEM began to grow in the engineering and physics communities, 
more applied mathematicians became interested in giving the method a firm mathematical 
foundation. Although the FEM has its inception in the late 1950’s as method for solving 
structural problems, it was not until 1963 that it got recognized as a general method that can 
deal with problems of the elastic continuum. Finally, Zienkiewicz and Cheung [66] recognized 
that the method is one of an even more general nature, applicable to all field problems which 
can be formulated as extremum or stationary value problems. 

Evaluation of the FEM has brought into existence a great body of literature including 
a text book by Zienkiewicz [67]. Although a major portion of the literature in the early 
1970’s deals with static and dynamic structural analysis, there has been a continuous steady 
increase in the number of applications in other fields. Unquestionably, the FEM is now a 
well-established and accepted engineering analysis tool. 

2.3 Variational Methods 

Often continuum problems have different but equivalent formulations-a differential equation 
and a Variational formulation. In the differential formulation, the problem is to integrate 
a differential equation or a system of differential equations subject to given boundary con- 
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ditions. The two problems are equivalent because the functions that satisfy the differential 
equations and their boundary conditions also extremize or make stationary the function- 
als. This equivalence is apparent form the calculation of variations, which shows that the 
functionals are extremized or made stationary only when one or more Euler equations and 
their boundary conditions are satisfied. And these equations are precisely the governing 
differential equations of the problem. The classical Variational formulation of a continuum 
problem often has several advantages over the differential formulation from the standpoint 
of obtaining an approximate solution. First, the functional, which may actually represent 
some physical quantity in the problem, contains lower order derivatives than the differen- 
tial operator, and consequently the approximate solution can be sought in a larger class of 
functions. Second, the problem may possess reciprocal Variational formulations, that is, one 
functional to be minimized and another functional of a different form to be maximized. In 
such cases we have a means for finding upper and lower bounds on the functional, and this 
capability may have significant engineering value. Third, the Variational formulation allows 
us to treat very complicated boundary conditions as natural boundary conditions. Fourth, 
form a purely mathematical standpoint the Variational formulation is helpful because with 
the calculus of variations it can sometimes be used to prove the existence of a solution. 

Historically, Variational methods are among the oldest means of obtaining solutions 
to problems in physics and engineering. One general method for obtaining approximate 
solutions to problems expressed in Variational form is known as the Rayleigh-Ritz method. 
Another one is the method of weighted residuals. 

2.3.1 Rayleigh-Ritz Method 

The Rayleigh-Ritz method consists of assuming the form of the unknown solution in terms 
of known functions (trial functions) with unknown adjustable parameters. From the family 
of trial functions we select the function which satisfies the functional theory. The procedure 
is to substitute the trial functions into the functional and thereby express the functional 
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in terms of the adjustable parameters. The function is then differentiated with respect to 
each parameter, and the resulting equations are set equal to zero. If there are ‘n’ unknown 
parameters, there will be ‘n’ simultaneous equations to be solved for these parameters. By 
this means the approximate solution is chosen from the family of assumed solutions. 

The procedure does nothing more than to give us the “best” solution from the family of 
assumed solutions. Clearly, then the accuracy of the approximate solution depends on the 
choice of trial functions. We require that the trial functions be defined over the whole solution 
domain and that they satisfy at least some and usually all of the boundary conditions. 
Consider, as an example the following simple one dimensional boundary value problem (bvp), 
consisting of the differential equation 

= fix) on 0<:r<X (2.1) 

ax ax 

u(0) = uo (2.2) 

Xp{X) = q (2.3) 

ax 

More generally, the differential equation is denoted by 

a(u) = / (2-4) 

Xet us denote the domain by Q and the boundary by dO,. The boundary condition given by 
equation (2.2) is usually called ‘Dirichlet boundary condition’. In the sequel, the part of the 
boundary where Dirichlet boundary has to be applied is denoted by 5f2o 

hiu) = go ( 2 - 5 ) 

The derivative boundary condition given by equation (2.3) is called ‘Neumann boundary 

condition’. In the sequel, it is denoted by dO^i- 

h{u) = 9i ( 2 . 6 ) 

The ’bvp given by equations (2.1) - (2.3) are said to be in its strong form, requiring the 
satisfaction of the differential equation (2.1) in all points of the domain dQ, the satisfaction 
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of the Dirichlet boundary condition given by equation (2.2) on Sflo the satisfaction 
of Neumann boundary condition equation given by equation (2.3) on 5Qi. In general the 
method is described below for a linear Variational problem. 

Consider the Variational problem of finding the solution u such that 

B{w,u) = l{vj) (2.7) 


for all sufficiently differentiable functions w satisfying the homogeneous form of any specified 
essential boundary condition on u. When the functional B is bilinear and symmetric and 
when I is linear, the problem in equation (2.7) is equivalent to minimization of the quadratic 
functional 

I{u) ='^B{u,u) -l{u) (2.8) 


In the Rayleigh-Ritz method, we select an approximate solution to equation (2.8) in the form 
of a finite series 


N 

'^N ~ ^ , ^j4*j d" ^0 (2.9) 

j=l 


where the constants cj, called the Ritz coefficients, are chosen such that the equation (2.7) 
holds for w = <l>i (i = 1,2,- • • ,N); i.e. equation (2.7) holds for N different choices of w, so 
that N independent algebraic equations in are obtained. The algebraic equation is 
obtained by substituting for w: 

= ^(^t) (* = 1, 2, • • • , Af) 

j=i 


If B is bilinear, the summation and the constants Cj can be taken outside the operator. We 
have 

N 

B{(i>^,6j)cj = - B{<f)^, 0o) or 

j=i 

N 

Y ByCj = Fi, Bij = B{<i)^, (pj), Fi = l{(t)i) - B{4>i, 4>o) (2.10) 


which represents the algebraic equation in a system of N linear algebraic equations in 
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N constants Cj. The columns (or rows) of the matrix By = must be linearly 

independent in order that the coefficient matrix given by equation (2.10) can be inverted. 
The requirements on and are as follows. 

In the approximation un in equation (2.9), we note that ujv must satisfy the specified essential 
boundary conditions of the problem; any specified natural boundary conditions are already 
included in the Variational problem in equation (2.7). The particular form of un in equation 
(2.9) facilitates satisfaction of specified boundary conditions. If we were to use the form 

TV 

Un = Cj(f>j{x) 
j=i 

then it would not be easy to satisfy non-homogeneous boundary conditions. For example, 
suppose that un is required to satisfy the conditions un{xo) = uq at a boundary point 

X = Xq: 

TV 

Cj^ji^o) =Uo 

- 3=1 

Since Cj are unknown parameters to be determined, it is easy to change (j)j (x) such that this 
relation holds. If uq = 0 then any such that 4>j{x) = 0 would meet the requirement. By 
writing the approximate solution un in the form of equation (2.9), a- sum of homogeneous 
and non-homogeneous parts, the non-homogeneous essential boundary conditions can be 
satisfied by (po, (poi^o) — uq, and (pj are required to satisfy the homogeneous form of 
the same boundary condition, (pjixo) = 0. In this way un satisfy the specified boundary 
conditions: 

TV 

WAr(a:o) = X) + (poixo) 

1=1 

= 0 - 1-140 

If all specified essential boundary conditions are homogeneous (i.e the specified value Uq 
is zero), then ^o{xo) is taken to be zero and <pj must still satisfy the same conditions, 
(pj = 0. Since <pi satisfy the homogenooiis essential boundary conditions, the choice w = <pi 
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is consistent with the requirement of a weight function. The approximate function (pi satisfy 
the following conditions, (i) (pi should be such that B{(pi,(pj) is well defined and non zero 
(ii) (pi must satisfy at least the homogeneous form of essential boundary conditions of the 
problem. The only role that (po plays is to satisfy the specified non-homogeneous essential 
boundary conditions of the problem. 

Further if the bilinear form is symmetric, the Rayleigh-Ritz method can be viewed as 
one that seeks a solution of the form given by equation (2.9) in which the parameters are 
determined by minimizing the quadratic functional corresponding to the symmetric bilinear 
form. If B is not symmetric, we do not have a quadratic functional. 

2.3.2 The Method of Weighted Residuals 

The first basic ingredient of the finite element method is that an approximate solution is 
sought which belongs to some finite dimension function space. In the weighted-residual 
method, the solution u is approximated, in much the same way as in the Rayleigh-Ritz 
method, by the expression: 

N 

u = 'tp+Y. (pkUk ( 2 . 11 ) 

k=l 

except that the requirement on ip and 4>k for the weighted-residual method are more strin- 
gent than those for the Rayleigh-Ritz method, ip is a, function which satisfies the boundary 
conditions given in equations (2.2) and (2.3) and ^k are functions that satisfy the boundary 
conditions of the same form but with right hand side equal to zero. For the given problem, 
the construction space (pk are generally called basic functions or shape functions. Since the 
dimension of the function space (p = {(pk : = 1, 2, • • • , iV} is finite and it yields a set of 

algebraic equations of the form: 


KU = F 


( 2 . 12 ) 
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where U is a vector consisting of elements {ui,U 2 , • • ■ ,ui), K is the coefficient matrix and F 
is a force vector. In general equation (2.11) can not satisfy the differential equation (2.1) at 
each point of the domain. This means that the approximate solution u can not be identical 
with the exact solution u. Of course, the shape functions should be chosen so that by 
enriching the function space i.e. letting N grow, the approximation obtained by equation 
(2.11) becomes better. This means that the approximate solution converges to the exact 
solution. This is called the completeness requirement of the function space. 

Since a function u given by equation (2.11) can not satisfy the differential equation 
(2.1), upon substitution of u into equation (2.4), a residual is left; 

rn = a{u) — f in f2 (2.13) 

An approximate solution of the bvp is now obtained by finding a way to make this resid- 
ual small in some sense. In the finite element method this is done by requiring that an 
approximate number of weighted integral of the residual over 12 be zero: 

[ wiradQ^O; 1=^1, 2, (2.14) 
Jn 

where W = {wi]l = 1,2,- ■■,N} is a set of weighting functions. Obviously, the conver- 
gence requirement now also implies a requirement of completeness of the space of weighting 
functions, i.e. equation (2.14) should imply rn — > 0 for N oo. 

Clearly with satisfaction of the completeness, for N oo, the weighted-residual for- 
mulation given by equation (2.14) for a function of the form equation (2.11) is completely 
equivalent to the strong formulation of the bvp defined by equations (2.1) - (2.3). An ap- 
proximate solution then is obtained for N being finite. 

2.3.3 Galerkin Formulation 

Among the possible choices for the set of weighting functions the following ones are the most 
obvious. The weighting functions can be chosen to be Dirac-Delta functions in JV points. 
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This choice boils down to making the residual equal to zero in a number of chosen points. 
This method is called the point-collation method. Obviously, it has much in common with 
the finite difference philosophy. 

A second possible choice of weighting functions is given by ' 


Wi 


1 for xi < X < xi+i 
0 for X < rc; or a; > xi+i 


The weighted residual statements, given in equation (2.14), now require the integral of the 
residual to be zero on N sub-domains. This method is called the sub-domain collocation 
method. 

The most popular choice for the weighting functions in the FEM is that the shape 
functions themselves are the weighting functions: wi = <pi. This method is called the 
Galerkin method. Its meaning is that the the residual is made to be orthogonal to the space 
of the shape functions. 

To illustrate the Galerkin method, consider the bvp (2.1), with constant A. Then: 

'Ip = Uo+^ X 


Consider, a Fourier-Sine expansion of u: 

u 


^ Tik'x 


Uk sin 

k=l 


X 


where k' = k — \ Then: 


^ irk' ^ . 'Kk'x 

rn = -X X. sin-Y- - f 


k=l 


The Galerkin method then gives 
^ ,'Kk' ^ 


A 53 ^fc(' 

fc=i 


A 


■) r 


-Kk'x . irl'x 
sin- -— sin-zT- ax 


X 


X 


-L 


sm-^ j ax 
A 


Then noting that: 



irk'x . ttZ'x 


s^n- 




sin- 


X 



for k' = I' 
for k' ^ r 
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we find: 


Ul 


2X 


/ ^ 


. ttI^x ^ 
sin--:— ax 

A 


In some situations, modification of the standard choice are necessary. The standard choice 

is denoted by Bubnov — Galerkin method. When modified weighting functions are used, the 

method is denoted by the term Petrov — Galerkin method, i.e. wi ^ (j)i. To illustrate the 

Bubnov — Galerkin and Petrov — Galerin methods, the following example is considered. 

(Pu o 


For a weighted residual method, 4*1 should satisfy the following boundary conditions 

(bcs). 

^o(O) = 0, 4>o'{l) = 1 (satisfy actual bcs) 

4>i{0) = 0, ^/(l) = 0 (satisfy homogeneous form of a specified bcs) 


For a choice of algebraic equations, we assume (f>o{x) = a + bx and use the two constants 
on <po to determine the constants a and b. We obtain (j)Q{x) = x. Since there are two 
homogeneous conditions, we must assume at least a three parameter polynomial to obtain a 
nonzero function <j)i = a + bx + cx"^. Using the constants on (pi, we obtain 

4)1 = — c x (2 — x) 

The constant c can be set equal to unity because it will be absorbed into any one of the 
parameters, we can assume one of the forms 4>2 = (^ + bx + dx^ ox 4>2 = o. + cx^ + dx^. With 
d ^ 0, 4>2 does not contain all order terms in either case, but the approximate solution is 
complete because {^i, 4>2} contains all terms up to degree three. For the first choice of (P2, 
we obtain 

2 

4>2 =x'^{l - -x) 
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The residual in the approximation of the equation is 

2 cP4,. 


R~ - (0 + ^'7^^ ~ (00 + X] ^*00 + ^ 

i=l z=l 

2 

= Ci(2 — 2x + x^) + C 2 (— 2 + Ax — x^ + -x^) — x + x^ 


In the Bubnov — Galerkin method, ‘ipi = (f>i, we have 

pi /‘i o 2 

/ x{2 — x) R dx = 0, / a;^(l — -x)R dx = 0 

Jo Jo 3 


or 


u = 1.2894X - 0.1392:2 - 0.003252;^ 


In the Petrov — Galerkin method, ^ <pi. Let the weighting functions be t/’i = x, "02 = 

we have 

[ xR dx = 0, [ x^R dx = 0 

Jo Jo 

or 

u = 1.302053a: - 0.173021a:2 - 0.014663a;^ 

2.3.4 Weak Formulation 

In many problems, it is not practical to construct a function which satisfies the boundary 
conditions in order to arrive at an expression for the approximate solution consisting of set of 
algebraic equation K U = F as discussed in section (2.3.2). More generally , an approximate 
solution can be expressed as 

N 

fc=i 

now the approximate solution not only has a residual with respect to the field equation (2.4), 
but also with respect to the boundary equations (2.5) and (2.6). 

ro = bo{u) - go 
n = bi{u) - gi 


(2.16) 

(2.17) 
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A weighted residual statement is now to be of the form 

I wiTq dCl + f WiTo ddO, + / wlri ddO, = 0 (2.18) 

J dQo J dVt\ 

This complicates the formulation since additional weighting functions on the boundaries are 
now to be chosen. Since the number of degrees of freedom of the approximate solution, 
equation (2.15) is N, an equal number of independent weighting functions Wi can be chosen, 
while and w] are to be dependent on wi. In practice very often, there is a natural way to 
chose the dependent weighting functions. For the bvp equations (2.1) - (2.3), above equation 
(2.18) becomes 

d dij di'i 

Jo ^^^dx ^^dx'^ -f]dx + w^[ u(0) - uq ] +wl [ A— (AT) - q ] (2.19) 

Of course, for this problem the weighting functions at the boundary reduce to weighting 
factors wf and wf. By integration by parts on the first term, this becomes 

■ i «-fdx + w,[{u(0)-uo] 

+ wj [ a£(X) - 5 ] = 0 (2.20) 

The weighted residual statement is simplified considerably by choosing the weighting factors 
on the Neumann boundary as 

w} = - wi{X) 

The weighted residual statement then becomes 

, dwi du , , du , , nr.. /..N 1 

- i ~ i “(»)-“« 1 

4- wi{X)q = 0 (2.21) 

Furthermore, if the Dirichlet boundary condition can be imposed on the approximate so- 
lution, the weighting functions and the weighting factors can be chosen to be zero at the 
Dirichlet boundary, so that the residual statement further simplifies to 

_ dx - wifdx + Wi(X) q = 0 (2.22) 

Jo ax dx Jo 
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subject to the Dirichlet boundary conditions 

u(0) = Uq 
wi(0) = 0 

The weighted residual statement in the form of equation (2.22) is called the weak formulation, 
Clearly, the weak formulation given in equation (2.22) along with the Dirichlet boundary 
conditions is not completely equivalent to the strong formulation even as N ^ oo. By the 
construction of the weak formulation, any solution of the strong formulation satisfies the 
weak formulation. The converse, however, is not true. The weak formulation allows so- 
lutions which have a lower degree of regularity than required for strong formulation. This 
is the origin of the terms weak and strong. For instance for the bvp equations (2.1) - (2.3), 
the solution must have a continuous first derivatives. We express this by saying that the 
degree of the regularity is to be C^. Whereas, the corresponding weak formulation, on- 
ly requires continuity of the function value itself; the necessary degree of regularity here 
is C^. This means that functions with discontinuous first derivatives are allowed in the 
weak formulation. 

Further, in the weak formulation the Neumann boundary condition need not be im- 
posed in an explicit way to the solution. Boundary conditions of this type enter through 
the integration by parts in a natural way into the formulation. Therefore, these boundary 
conditions are called natural boundary conditions. Also, the boundary conditions which 
have to be imposed explicitly in the weak formulation are called as essential boundary 
conditions. 

2.3.5 Shape Functions 

A second basic ingredient of the finite element method is the piecewise manner in which 
the shape and weighting functions are constructed. The domain Vt is subdivided into non- 
overlapping subdomains fig, called “elements simplified geometrical form”. In general, the 
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integral in the weak f ormulation can be split into a sum of integrals over elements. 

/ (F) dCl ^ Y. L (F) dCl 

Then, obviously, in the piecewise contributions to the integrals, it is computationally advan- 
tageous to have as many zero contributions as possible. This is achieved when the shape 
functions and the weighting functions associated to some index are only non-zero in the 
smallest possible number of the elements associated to this index. Shape and weighting 
functions which are only non-zero in a small set of elements are said to have a compact 
support. In the finite element method, shape and weighting functions with compact support 
are constructed from an interpolation problem over the domain. For instance, a function u 
which is obtained through linear interpolation between function values defined in the grid 
points of the grid can be written as 

L 

^=13 4>kUk (2.23) 

k=0 

The shape functions (fik in this expression look like a hat form. For a function representation 
based on an interpolation, the values Uk in the equation (2.23) have the meaning of function 
values in the gird points. Obviously, other interpolation schemes are possible. For instance, 
the function u could be obtained by piecewise constant interpolation or the interpolation 
could be piecewise quadratic. 

In all these cases, the values of Uk represent function values in some points associated 
to the elements. In the finite element technique, these points are called nodes. As a simple 
example, we first consider the bvp (2.1) - (2.3) with piecewise linear shape functions and a 
standard Galerkin weak formulation. The approximate solution is then represented by 

* = Z] ^kUk = Y. (2-24) 

k=0 e j=l 

In equation (2.24) the sum over the nodes is rearranged as a double sum, first over the 
elements and then over the nodes belonging to the element. The shape function, (pk, associ- 
ated with the nodes, are called the global shape functions. On the element level, the shape 
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functions are called local shape functions or element shape functions. The linear shape 
functions are shown in Figure 2.2: (a). 



(a) 


(-1. 1) (1, 1) 


T 



*"4 


(- 1 ,- 1 ) 


(b) 



(c) 

Figure 2.2; (a) Linear one dimensional elements (b) Bilinear rectangular elements and (c) Its shape 
functions in normalized coordinates. 


For the linear shape functions: 


Hence: 


^ AXe 

II 

d4>l _ -1 

d(f)2 

dx AXe ’ 

dx 


X — Xe- 
AXe 


1 

AXe 


1 


The Galerkin weak formulation given in equation (2.22) is 

[ ^ dx = - [ wifdx + wi(X)q 

Jo dx dx Jo 


(2.25) 
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Although by the piecewise shape representation the approximate solution can not satisfy the 

strong form of the boundary value problem in any point of the domain, the solution obtained 

from the weakformulation represents a valid approximation of the problem. 

A rectangular bilinear element is a four noded serendipity class of element. A typical 

bilinear element in local coordinates is shown in Figure 2 . 2 : (b). The linear variables 

in ^ and 77 directions along the edges are defined by ^0 = ^ 6 = 1 or — 1 and 770 = 

V Vii ??i = 1 or — 1 . The linear nodal shape functions in ^ and 77 directions at the nodal 

points ^=lor — 1 , 77 = 1 or — 1 are defined by 

(l + '^^o) jl+VVo) 

2 2 


respectively. Now for the bilinear element the nodal shape functions can be derived by taking 
the product of the two linear shape functions associated with that node. The bilinear nodal 

shape function associated with node i (Shown in Figure 2.2: (c)) are given by 

^ _ (l + -^?o)(l + Wo) o ^ . 

cpi — - j % — 1 5 2 , o 5 4 


The element is termed bilinear because of the product of ^ and 77 in addition to the linear 
terms in the expression. The form given for 4 >i automatically satisfies the usual requirement, 
i.e. the value at node i is unity and at all other nodes of the element is zero. 


2.4 Upwind Finite Elements 

Initial finite element formulations for convective transport problems also used the Galerkin 
method, but with mixed results. In fluid flows or convective heat transfer, the matrix asso- 
ciated with the convective term is non symmetric, and as a result, the ‘best approximation’ 
property is lost. In practice, solutions are often corrupted by spurious node-to-node oscil- 
lations or ‘wiggles’. Wiggles are most likely to appear in convection dominated cases (high 
Peclet or Reynolds number) when a downstream boundary condition forces a rapid change 
in the solution. An example of this phenomenon is shown in the Figure 2.3 for fluid past a 
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Figure 2.3: Oscillations generated upstream of a block for Reynolds number = 200 

block in a narrow channel [68] . The only way to eliminate the oscillations is to severely re- 
fine the mesh such that convection no longer dominates on element levels. Obviously, mesh 
refinement is needed in regions where accurate representation of the solution is required. 
However, often only the global solution features are desired and in this case mesh refine- 
ment is required only to prevent oscillations. This provides motivation for development of 
an alternative to the Galerkin formulation which precludes spurious oscillations, regardless 
of mesh refinement. It is also well known that the Galerkin finite element method give rise 
to central-difference type approximations of diflferential operators. It is thus not surprising 
that ‘wiggles’ have also affected central-difference solutions. In fact, obtaining solutions at 
all has been problematic due to the reliance on iterative solvers in finite difference coupled- 
difference equations. It was discovered, however, that wiggle free solutions could be obtained 
by the use of ‘upwind’ differencing on the convective term. Upwind differencing amounts 
to approximating the convective derivative with solution values at the upstream and central 
nodes of a three node stencil. It is well known that the upwind convective term can be con- 
structed simply by adding artificial diffusion to a central difference treatment. This artificial 
diffusion interpretation has been the basis for extensive criticism of upwind methods. In a fi- 
nite element framework, upwind convective terms can be developed in several different ways. 
The initial upwind finite element formulation, for the one dimensional advection-diffusion 
equation, employed modified weighting functions to achieve the upwind effect. In essence, 
the element upstream of a node is weighted more heavily than the element downstream of a 
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node. 

Many of the optimal upwind finite element formulations lead to the same system of 
matrix equations and give exact solutions for the one-dimensional model problem. Unfor- 
tunately, when generated to more complicated situations, some of these formulations are 
far from optimal. In multidimensional problems, solutions often exhibit excessive diffusion 
perpendicular to the flow direction. Overly, diffuse results also appear in transient problems, 
or when source terms are present. In addition, in many instances the Galerkin formulation 
gives oscillation-free solutions which are more accurate than upwind solutions, obviously, 
such results cause upwind techniques to be viewed with some suspicion (and even contempt). 
Some investigators believe that use of the Galerkin formulation, along with mesh refinement, 
is the only possible way to achieve solutions that are both accurate and wiggle free. 

To illustrate upwind finite elements, consider the following example. Let Q be a bounded 
region in Ji", n > 2, with piecewise smooth boundary. A general point in Q, is denoted by 
X. The model equation under consideration is the steady advection-diffusion equation in an 
incompressible flow field: 

^74) = V^kVcj) (2.26) 


where u is the flow velocity, k is the diffusivity and V is the gradient operator. We assume k 

to be positive definite (both u and k are given functions of x). The boundary value problem 

for equation (2.26) consists of finding a function (j), satisfying equation (2.26), such that: 

4> = f on Ti 1 

kV 4) = g on La j 


(2.27) 


where / and g are given functions of x, n is the outward unit normal vector to F, and F i and 
F 2 are subregions of F satisfying Fi UF 2 = F and Fi n F 2 = 0. A weak form of the problem 
is given by 

[ {wu'^V4> + Vw'^kV4> ) dQ = [ w gdn (2.28) 

Jfl ‘'^2 


where 4) is required to satisfy the first of equation (2.27) and w is required to be zero on Fi. 
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The second of equation (2.27) is a natural boundary condition implied by equation (2.28). 
The one dimensional version of equation (2.26) is 


d4> ___ d ^d(f) 
dx dx dx 


(2.29) 


where u and k are given functions of x, k is positive, and fi = (0, L). Various specializations 


of boundary conditions are possible. For example, 

m=fo \ 
m = fL I 


(2.30) 


or 


m = A 
m i(i ) = gi 


(2.31) 


are allowable combinations, where /o,/i, and gi are given constants. Weak form of the 
one-dimensional problem, corresponding to equations (2.29) - (2.31), are given by 




dd) dw , dS , , 

w u — + — k — ) dx = 

cLoc cLoc dx 


0 

w{L) gi 


(2.32) 


in the former case, 0 is required to satisfy equation (2.30) and rt;(0) = w{L) = 0, whereas 
in the latter case, (f) is required to satisfy only the first of equation (2.30) and u;(0) = 0. 
Approximate solution of equations (2.28) and (2.32) may be constructed by the finite element 
method in which 

(f> ^ Ni di (2.33) 

i 

where is the shape function associated with node H’ and di is the approximate value of (j) 
at node H\ In the standard- Galerkin method (Bubnov- Galerkin method), w is taken to be, 
in turn, equal to each shape function which vanishes on Fi. In the one-dimensional case in 
which u and k are constant and element lengths are equal, employing piecewise linear shape 
functions, leads to a central difference approximation to the advection-diffusion term. Under 
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these circumstances it is well known 


that unless 
2k 


h < 


u 


(2.34) 


where h is the element length, the discrete solution will contain spurious oscillation called 
wiggles (move form side to side with rapid short moments). Analogous conditions to e- 
quation (2.34) holds in more general sense in which the Galerkin finite element method is 
employed and result in excessive mesh refinement in problems of interest. 

To circumvent such difficulties upwind finite element formulation have been developed by 
a more general weighted residual method is which w is chosen from a wider class of functions 
than the shape functions. The procedure, although effective leads to grater evaluations than 
those used in the Galerkin formulation. 

In general, there has been three basic techniques utilized to achieve the upwind effect 
in finite elements. 


1. ARTIFICIAL DIFFUSION. Artificial diffusion, is added to the physical diffusion, and 
a conventional Galerkin finite element discretization is employed. This is actually a 
‘balance diffusion’, in that it balances the negative diffusion of the Galerkin treatment. 

2. QUADRATURE. The numerical quadrature rule for the convective term is modified 
to achieve the upwind effect. 

3. PETROV-GALERKIN. The weighting functions for a typical node is modified to 
weight the element upwind of the node more heavily than the download element. An 
example is shown in Figure 2.4. 


2.5 SUPG strategy 

The SUPGFEM has the robust qualities of a classical upwind method, but is not subject to 
any of the artificial diffusion criticisms. The basic idea of the Streamline Upwind method is 
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GALERKIN 


FLOW 



Figure 2.4; Galerkin weighting function and Upwind Petrov- Galer kin weighting function. 

to add diffusion (or viscosity) which acts only in the flow direction. Extended to a Petrov- 
Galerkin formulation, the standard Galerkin weighting functions are modified by adding 
streamline upwind perturbation, which again acts only in the flow direction. The modified 
weighting function is applied to all terms in the equation, resulting in a consistent weighted 
residual formulation. 

Let be a bounded region in B? and assume has a piecewise smooth boundary P. 
Let X = {rri}, i = 1,2. denote a general point in SI and let n = {nj be the outward normal 
vector to T. Let Tp and be subsets of P which satisfy the following conditions: 

r7^ = r 

P3nP/i = ?i' 

The superimposed bar in the above equation represents set closure and (j) denotes the empty 
set. Consider a discretization of SI into element subdomains SI®, e = 1,2.. Let P® denotes 
the boundary of f2®. Further, assume 

Ue 0® = SI, He S^® = (^ 
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Finally, define the ‘interior boundary’, Fj^j as follows: 

= Ue - r 

For the illustration of the strategy, consider the unsteady linear advection-diff’usion equation 
for incompressible flow fields. 

t + cTi^i = f (2.35) 

where 

cTj = cr“ + af (total flux) 
cr“ = Ui(f> (advection flux) 

= — k^J 4>^ j (diffusive flux) 

In the above, / is a source term, Ui is the flow velocity, and % is the diffusivity. Each of /, 
Uj, and kij are assumed to be given functions of x and t. Further, the velocity field, Ui , is 
assumed to be divergence free. The initial- value problem defined by equation (2.35) consists 
of finding a function (?i(x, t) which satisfies equation (2.35) on Q, and 

4> = g on Fj 
- a^ = h on Fa 

^(X, 0) = (j)Q 

where g and h are given functions of x and t, and the initial condition is a given function 
of X. In a usual Galerkin weighted residual method, the weighting functions are considered 
to be continuous across inter element boundaries. The Streamline Upwind/Petrov-Galerkin 
formulation, however needs discontinuous weighting functions of the form 

W = w + P (2.36) 


where tu is a continuous weighting function, and P is the discontinuous upwind contribution. 
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Both w and P are assumed to be smooth on the element interiors. Consider a point x in 
V^nt■ Designate (arbitrarily) one side of Tint to be the ‘plus side’ and the other to be the 
‘minus side’. Let n'^ and n~ be unit normal vectors to Tmt at x which point in ‘plus’ and 
’minus’ directions, respectively. Clearly = n+. Let Tf and denote the values of Dj, 
obtained by approaching x from the positive and negative sides, respectively. The ‘jump’ in 
r„ at X is defined to be 

[ (7n ] = ( (Tf - a- )nf = atnf + 

it is also to be noted that ‘jump’ is invariant with respect to reversing the ‘plus’ and ‘minus’ 
designations. We shall also assume that the trial solutions, (p satisfy p = g on. Tg and 
weighting functions, w, satisfy tt; = 0 on F^. 

The Streamline Upwind /Petrov-Galerkin weighted residual formulation for the initial 
boundary value problem defined by equation (2.35) is 

[ w {(p + ) dO. - f w^icf dQ + [ P {p + - f) dCl 

Jn Jo. g Jo^ 

= [ w f dO, + [ w h dr. 

Jo Jr>^ 

Integration by parts yields 

Y f W {<p + ai^i - f) dn ~ [ w {ai + h) dV 

e J0‘ JTh 

- [w [af] dV = 0. 

As shown in the Figure 2.4, upwind finite elements may be constructed via the Galerkin 
method with added artificial diffusion. In the one-dimensional case, the artificial diffusion, 
when optimally selected, balanced the negative diffusion inherent in the Galerkin method, 
resulting in exact nodal solutions for the model problem. However, the multidimensional 
generalizations like flow of lubricant in parallel slider bearings of some finite element upwind 
schemes were unsuccessful due to the cross wind diffusion problem. It was apparent that the 
upwind effect, arrived at by whatever means, was needed only in the direction of flow. In 
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As shown in the equation (2.36) the perturbation is seen to be 

p_ kUj Wj 

M 

Clearly, if w is continuous across element boundaries, then P, and thus W, will be dis- 
continuous. Examples of w and W for one-dimensional linear finite elements are shown in 
Figure 2.5. 


Flow 



Figure 2.5: Comparison of Streamline Upwind/Petrov-Galerkin (SUPG) and Galerkin weighting 
functions for a node A 


2.6 Implementation of the Finite Element Method 


The final ingredient of the finite element method is the way in which the nodal equations 
are constructed. As seen in the example, the integration involved in the weak formulation 
leads to a sum of contributions on elements adjacent to the node which is treated. Instead 
of performing the integration on the set of elements adjacent to a node we could of course 
perform the integration on all elements separately and then afterwards construct the nodal 
equations by adding contributions from adjacent elements. This is called assembly. It has 
the advantage that the integration on all elements can be done in the same way and, as a 
consequence, with the same routine. For instance, for the equation (2.22), the integral in 
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An element matrix can be defined by 

= A / 




dx dx 


dx 


The component of the system matrix in the global discrete system 

KU = F 


( 2 . 37 ) 


can now be found by adding components from the element matrices. The system matrix K 
is called stiffness matrix. In general problems, the computation of the element stiffness 
matrix and the element right hand side can not be performed analytically due to complexity 
of the integrand and due to the complexity of the element. In such cases it is generally 
advisable to go directly for an iterative technique. 



Chapter 3 

Lubrication in Tilted Pad Slider 
Bearings: Effect of Flux Boundary 
Conditions on the Pad^ 


3.1 Introduction 


A tilted pad slider bearing consists of two surfaces (as shown in the Figure 3.1 separated by 
a lubricant film. One surface is stationary while the other moves with a uniform velocity. 
The direction of motion and the inchnation of planes are such that a converging oil film is 
formed between the two surfaces and the physical wedge pressure generating mechanism is 
developed in the lubricant film. It is this pressure generating mechanism that enables the 
bearing to support a load. Several researchers have considered this basic geometry under 
various situations. 

However, the parallel sliding has not been investigated so thoroughly. Parallel sliding 
of surfaces in the presence of a lubricant represents a common bearing configuration which 
can be found, in many engineering applications like Mechanical seals. Machine tool ways. 
Piston rings. Plain collar thrust bearings etc. 

Several researchers like Fogg [36] , Osterle et al. [38], Lewicki [39], Cameron [43] , 

^Accepted for publication in the “International Journal of Numerical Heat Transfer” 
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Young [69], Lebeck [47, 48] etc., have experimentally investigated the thermal influence on 
load carrying capacity of parallel slider bearings. Recently Rodkiewicz and Sinha [58] and 
Schumack [19] have numerically analyzed the problem of thermohydrodynamic lubrication 
with Dirichlet boundary prescription on the pad and slider. But the actual slider bearing 
configuration would require zero/nonzero flux boundary prescriptions. Gero and Ettles [18] 
has attempted the problem but his study was focussed on scheme development. Further he 
has dealt with either constant viscosity /density or linearized viscosity with constant density. 
The influence of thermal boosting, which amounts to considering slider, pad and inflow 
lubricant at different temperatures, on load carrying capacity, drag force etc of the lubricant 
in thermohydrodynamic slider bearing context has not been investigated. 

So the aim of the study in this Chapter is to numerically investigate the influence of 
various temperature boundary prescriptions and thermal boosting on velocity distribution, 
temperature distribution, pressure distribution, load carrying capacity and frictional drag 
force in the lubricant film by analyzing the non-linear governing equations with tempera- 
ture dependent viscosity and density at industrially significant (Rodkiewicz and Sinha [58], 
Gethin [17]) flow speeds and inlet pressures. For the numerical simulation of hydrodynam- 
ic lubrication problems Reddi [70] has proposed a variational principle based on Reynolds 
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Equation. Later Huebner [14] has shown the applicability of finite element method for 
thermohydrodynamic lubrication analysis based on variational principle for pressures, and 
Bubnov-Galerkin weighted residual approach for energy equation. He assumed the viscosity 
and density to be temperature dependent and obtained results considering high inlet pres- 
sures and high speeds. The scheme proposed by Huebner [14] suffer from well known node 
to node spurious oscillations. As we discussed in Chapter 2, Upwind strategies are known to 
preclude such oscillations in the solution. Brooks and Hughes [68] proposed SUPG Galerkin 
formulation for convection dominated flows with particular emphasis on the incompressible 
Navier-Stokes equations. SUPG Finite Element Method has been successfully used for ther- 
mal flow analysis by Grygiel and Tanguy [71] in advection diffusion problems. In the present 
Chapter the problem of thermohydrodynamic lubrication of slider bearing with temperature 
dependent density and viscosity which brings in extra pseudo compressibility terms into 
the energy equation has been successfully solved by using SUPGFEM. The results from the 
numerical simulations have been presented in computer generated plots and in Tables and 
the method has been validated for a hydrodynamic lubrication where analytical solution is 
available. 

3.2 Governing Equations 

The geometry and the coordinate system for a slider bearing analyzed in this Chapter is 
shown in Figure 3.1 Note that the dimensions are not shown to scale, the length of the 
bearing is typically of the 1000 times larger than the distance between the two surfaces. The 
basic assumptions for the lubrication flows are as follows. 

• Pressure is invariant across the fluid film (in the y direction) 

• Inertia forces are negligible 

• Velocity gradients in all but the y direction are negligible 

• The height of the fluid film y is very small compared to the span and length x. This 
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permits to ignore the curvature of the fluid film 

• The flow is laminar., no vortex flow and no turbulence occur any where in the film. 

• No external forces act on the film 

• No slip at the bearing surfaces 


In addition, the flow is assumed to be steady and laminar The Navier-Stokes equations thus 


reduce to 


9p_d du 
dx dy 


(3.1) 


and the continuity equation is 


9 ^ 


dy 


ifyu) = 0 


(3.2) 


where, p = p(T), p, = p{T). The boundary conditions are: 

u = U and u = 0 at y = 0, u = u = 0 at y = h (3.3) 


Combining equations (3.1), (3.2) and applying the boundary conditions given in equation 

(3.3) leads to the steady w . Reynolds equation for the pressure distribution. 


d dp. ^dh 
dx 6pdx dx 


(3.4) 


The assumptions for the energy equation are as follows. 

• Conduction terms other than those across the fluid film (in the y direction) are negli- 
gible. 

• Thermal conductivity and specific heat are constant. 

These assumptions combined with previous assumptions lead to the following form of the 
energy equation. 


^ dT dT, , d^T , dp 


dx 


(3.5) 


Viscosity and density are related to temperature via the following relationship. 

p = [1 - \(T - Ta)|, = P. exp[-^(r - TJj 


(3.6) 
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Now the following non dimensional variables (Schumack [19]) are introduced. 


u 


+ 


n 


y+ = VB_ 

^ Uho' 


““ B’ 


y+ = JL 


/lo’ 




Ta 


(subscript “a” indicates ambient condition), 


/i- = 


— JL ri+ — _£. 


/ia ’ 


Pa^ 


p'^ = (p is the gage pressure), 


_ _pJiL 

flaUB 


o B ’ 

A+ = XT,, PrE, = 


/i+ = %, /i+ = = /3T„ 


ho^ 

P — Po^/c/l^ 
— jbB • 


With this non dimensionalization, the governing equations are rewritten as follows. 

d dp-^ dh-^ 

Pi \-y “ 


dx^ ^ 6/i+ dx~^ ^ dx~^ 


(3.7) 


with the boundary conditions 


^Pt = 0 and = 0 at x 


— 


.+ — 




(3.8) 


and the corresponding energy equation in non dimensional form is as follows: 

Id^T+^PrE, + du+..PrE, +dp+ 


/)+ (u+ 




■u' 




(3.9) 


where p'^ and p'^ are given by: 

= 1 - A+(T'^ - 1), = exp[-P+(r+ -1)] 


(3.10) 


In view of the various industrial applications, the energy equation has been solved under 
the following different temperature boundary conditions: 

• Casel: Full Dirichlet model 
T'‘'(x+, 0) = T+, on slider, 

T+(x+, h+) = T+, on pad, where, = T/ = ^ 

• Case2: Zero half flux model 
T+( 2 ;+, 0) = T+, on slider, 

^{x*.h+) = 0, on pad (where n is the unit outward normal vector). 

• Case3: Non zero half flux model 
T+(x+, 0) = T+, on slider. 
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^ 0, on pad. 

• Cased: Zero full flux model 
%-(a:+,0) = 0, on slider, 

/i+) = 0, on pad 

• Case5: Thermal boosting model 
T+(x+,0) = r+, on slider, 

T'^{x+,h+) = Tp+, on pad with, 

(a)T+ = T+>i;+ (b) rp+ > T+ > 1^+ (c)r+>r+>j;+ 

(d) T+ = r+ < Tt (e) r+ > r+ < r+ (f) t+ > t+ < 

Actually, the combination of T^, T^, which enhances the load carrying capacity of a slider 
bearing constitutes a thermal boosting model. At the inlet/outlet the following boundary 
conditions are used. 

y^) = at the inlet boundary and (1, yA) = 0 at the outlet boundary. 

The load carrying capacity {W'^) and the frictional drag (F"'') have been determined from 
the following expressions. 

W7+ = 

F+ = 

where, W+ = and F+ = ^ 

3-3 Numerical Formulation of the Problem 

3.3.1 Variational/ Weak Formulation 

Let Q the domain of interest be bounded in with a piecewise smooth boundary T . Let 
(a;+, y+) denote the vector of spatial coordinate of a generic point in 0 and n is the outward 
unit normal vector drawn at any point (a;+, y+) on T. Now the variational principle n(p) for 
the Reynolds equation (3.7) in Q is derived as follows. Define the first variation in the integral 


f dx^ 

i: 


du'^ , 

y +=0 


dx'^ 


(3.11) 

(3.12) 




’txr 




IPU 
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form: 




(3.13) 


And the stationary of 11 (p) with respect to variations in p among the set of functions 
satisfying the boundary conditions requires that 


Jn 




dQ 


~ jjP 




Simplifying further 


dx^ ^ dx -^ ' 

ft+® dp+ a r d r d 


(3,14) 


5a:''' 


dp'*' 5a;''' ^ io dp'*' 5a;+ 50;-^ 




Now using the Green-Gauss identity and boundary conditions, the above equation leads to 
the following variational principle. 






dx'^ 


dQ, 


(3.15) 


where, p"^ and p"*" are defined in equation (3.10). The Bubnov-Galerkin weighted resid- 
ual formulation of convection dominated flows are known to lead to spurious node to node 
oscillations in the solution. A similar situation has been encountered with the energy equa- 


tion in the present study too. In this study to preclude such oscillations in the solution, 
the weighted residual formulation of the energy equation based on SUPG weight function is 
considered. The modified weighted residual formulation of the energy equation is given by: 


Wi[p+(u+ 


dT+ ld^T+ PrEc^,du+,^ PrEc__^dp+ 


dx^ 


+ V 


dy^ ^ Pe 




r-^c 4 - 


dx"^ 


; dn = Q (3.16) 


where, Wi is the SUPG weight function. It is to be noted that p+, p+ are dependent on 
T'^ and are given in equation (3.10) 


3.3.2 Finite Element Formulation 

The finite element discretization procedure reduces the lubrication problem to one of finite 
number of unknowns. This is accomplished by dividing the solution domain into bilinear 
rectangular elements and by expressing the field variables (pressure, velocities and tern- 



S.3 Numerical Formulation of the Problem 


59 


perature) in terms of approximating trial functions (Ni) within each element. Now is 
discretized into Ne bilinear rectangular elements such that: 

= n Q^ = (i) (3.17) 


where, denotes the interior domain of an element. Let be the boundary of 0®. The 
discretized representation of the field variables are given by: 

^el ^ N^l N^i TV’s/ 

^^kuf , < ~ E , pt ~ E ^kT^° • (3-18) 

^=1 k=l h=l jfc=l 

where, N^i = 2 for pressure case and Nei = 4 for velocities and temperatures. On 
introducing these elemental discretization details into equation (3.15) and (3.16), we would 
obtain element level representations. 

The element level representation of equation (3.15) is given by: 



On minimizing the Reynolds equation (3.15) with respect to variables pj, That is. 

duip) 


dpt 


0. 


dUjp) 

dpi 


dpi 




h+^ , d 


d 




5l2/i+ dx-^ ^ 


1=1 


On further simplification, we arrive at the following equation: 

/• h+^ FIATe r Qff? 


■* 'et /* 

si' 


h+^ dN^dNf. 
6/1+ 5a;+ 5a:+' 


_Le 

Pt 


dQ.^ 


= [ 


h+ 


dx+ 




(3.20) 


1, 2, 3, ■ • • , Ne/- On expanding the above equation and Nei — 2 for nodal pressures, we 


get the following elemental matrix equation. 

h+ dNxdNx woe f dElMl WQe 

in' 6u+ dx+ ax+ ““ -in® 6 a4+ dx+ dx+ 


In 

li 


n« 


6^“^ dx^ 

h+ MzMl woe 
6 m + dx^dx+ 


// 




h+ woe 

6fi+ dx+ dx+ ““ 


pi 

p2 


/ne h+p dQ' 

In^ 


(3.21) 


To evaluate the integrals in the above matrix equation (3.21), we have used, Simpson’s 
rule. Thus, solving the system of matrix algebraic equations using the well known 
Gaussian elimination scheme we obtain the nodal pressures. In the iteration process for the 
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first iteration is a constant. But in the subsequent iteration onwards it is a function of 
temperature defined by the relation given in the equation (3.10). Since the bearing solids are 
assumed to be rigid and the lubricant thermal effects are to be incorporated, it is essential 
that the bearing analyst must take into account viscous heating in the lubricant film and 
the resulting changes in lubricant viscosity. Such kind of analysis must include not only the 
generalized Reynolds equation, but also the conservation equation of thermal energy. The 
Reynolds and the energy equation must be solved simultaneously using SUPG finite element 
method. 

Now with the SUPG analysis of Brooks and Hughes [68] and Grygiel and Tanguy [71] 
one obtains the following simplified discretized elemental energy equation: 

E j {W{l-XU'£NtTr-l)) [(ExtX*) 


dNf 


dx'^ 


^ + (E<m) ^ ] 


dy+ 




dy+ dy+ 


^(exp(-^+(i:jv|Tr - 1))) 

k k 

(E {Acz mpf))] (3.2 

k k 


With W = Wi + Pi and Wi = Ni Pi = The upwind parameter k is calculated 

using elemental dimensions and elemental velocity (Grygiel and Tanguy [71], Brooks and 
Hughes [68]). 


3.3.3 SUPG Algorithm for Energy Equation 

Brooks and Hughes [68] SUPG algorithm has been incorporated for the present investigation. 

• stepl: Calculate (u+,-u+) at the element center. 

• step2: Calculate e^, e,, (unit vectors) 


• step3: Calculate Vj, where = e^.u, Vr, — e^.v 

• step4: Calculate h^, hr, (elemental characteristic lengths) 

• steps : Calculate a^, Or, where cxn — 
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• step 6 : Calculate fj where | = (Coth a{) — ^ and fj = (Coth ^ 

• step7: Calculate k where k = 

• step 8 : Calculate Pi where Pj = ^ m (sum). 

Now based on equations (3.21), (3.22) the numerical simulations have been carried out on a 
30 X 30 finite element mesh as shown in the Figure 3.2. The results have been obtained to 
an accuracy of e = 5.0 x 10“^ i.e. 

I ^New _ ^Old I 

)> 2 ' 

for j = l,2 ,---,n,--- 

where 7 ^ represents any of the field variables T/. 



Figure 3.2: Domain discretization for 30 x 30 mesh system 
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Figure 3.3: Flow chart for solution methodology 
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3.4 Results and Discussion 

The bearing characteristics in the present study appear as functions of A"*", k and dimen- 
sionless parameters Pr, Pg, Ec ■ For the present investigation, the value of these parameters 
are chosen in accordance with Lebeck [47, 48] and are as follows. 

U = 20 m/s, hi = 0.00005 m, B = 0.1 m, c = 1926 j/kg K, T„ = 310 K, pa = 
897.1 kg/m^, Pa = 0.0174 Pa.s, /3 = 0.035/ir, A = 0.0012/iif. The simulations have 
been carried out for various values of k,pf, and with the fixed values for the parameters 
Pj. (Prandtl number), Pg (Eckert number) and Pg (Peclet number). The results have been 
analyzed for load carrying capacity, drag force and velocity, pressure, temperature fields 
and are presented in the form of computer generated plots and Tables. To ensure the grid 
independence of the results, the numerical simulations were carried out on different mesh 
systems consisting of 10 x 10, 20 x 20, 30 x 30 and 40 x 40 bilinear rectangular elements. The 
results for load carrying capacity for different mesh systems are shown in Figure 3.4. 



Figure 3.4: Load carrying capacity on different mesh systems 

From the Figure 3.4, it can be easily concluded that the 30 x 30 mesh system yields 
a grid independent solution. Further the results have also been validated by comparing the 
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load carrying capacity with the analytical solution for the hydrodynamic lubrication. These 
comparisons are given in the Figure 3.5 for various inclinations of the slider bearing, i.e. 
with k = 0.4, 0.6, 0.8, 1.0 and at pf =■ 1, 10, 50 . Here “analy” stands for analytical solution 
and “nume” stands for numerical solution. 



Figure 3.5: Comparison of analytical and numerical load carrying capacity 

The dimensionless pressures at different locations in the sliding direction with k = 0.4 
at different inlet pressures have been shown in Figure 3.6. Also for the cases k = 0.8, k = l 
the non dimensional pressures in the direction of the sliding for values oipf = 1, 2, 5, 8 have 
been plotted and are shown in Figures 3.7 and 3.8. 

It is clearly seen from the Figures 3.6, 3.7 and 3.8 that for a fixed value of a;+, as the 
inlet pressure is raised the fluid pressure becomes higher. As a consequence the load carrying 
capacity is expected to be higher when the inlet pressure is higher for all values of k. 

From Table 3.1, it can be seen that consideration of non-zero inlet pressure leads to 
increased load carrying capacity. Further with Reynolds equation for pressure and temper- 
ature dependent viscosity, different flux boundary conditions do not alter the load carrying 
capacity of slider bearings with fixed values of k (inclination of the pad) and pf (inlet 
pressure) significantly. 




Nondimensional pressures, p* Nondimensional pressures. 



x+ 


Figure 3.7: Pressure distribution for the case k = 0.8 
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Figure 3.8: Pressure distribution for the case k^l 


Table 3.1: Coruparison of load carrying capacity with k for different flux boundary conditions 
at: i;+ = T+ = T+ = 1.5 


k 

0.4 

0.6 

0.8 

1.0 

with = 0 

PF'''‘(no flux) 
T^+(half flux) 
W^‘''(full flux) 

0.0003886 

0.0003887 

0.0003881 

0.00016031 

0.00016041 

0.00016031 

0.00005461 

0.00005428 

0.00005405 

0.000005785 

0.000005784 

0.000005651 

with pI = 1 
flux), 
W+(half flxix) 
^^■*■(^11 flux) 

0.7108 

0.7104 

0.71047 

0.6224 

0.6239 

0.623939 

0.55424 

0.554245 

0.5542461 

0.5000002 

0.5000002 

0.5000002 

with pI = b 
W'*'(no flux) 
W+(half flux) 
^’■^(full flux) 

3.5387 

3.5385 

3.53857 

3.1056 

3.1055 

3.10554 

2.76905 

2.7689 

2.76898 

2.4999 

2.4999 

2.4999 
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In Figure 3.9 pressure generated along the pad for k = 0.6, 0.4 at inlet pressure 
p+ = 0 have been plotted. 



Figure 3.9: Pad pressures at pf = 0 


In both the cases the peak pad pressures are noticed closer to the outlet. The quantita- 
tive behavior is in conformity with the existing results for a slider bearing. It is also further 
observed from these Figures 3.6 - 3.8 that the pressure profiles are nonlinear for k = 0.4 
and k = 0.8, whereas for k = 1, they are almost linear, indicating that no appreciable fluid 
pressure is generated for parallel slider bearing. Prom Table 3.2, It is observed that, lubri- 


Table 3.2: Effect of /5+ on load carrying capacity at different inlet pressures 



10.0 

8.0 

6.0 

4.0 

2.0 

0.0 

Pi =0 

0.0003886 

0.0011669 

0.0035013 

0.1049528 

0.31444689 

0.95129 

pt = 1 

w+ 

0.71086720 

0.718673 

0.7420864 

0.8123493 

1.022427 

1.658269 

pf = 5 

w+ 

3.5387 

3.546638 

3.570199 

3.641908 

3.85248 

4.486161 
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cants with low values of viscosity coefficient {^'^) enhance the load carrying capacity of a 
slider bearing irrespective of the inlet pressure. According to equation (3.10), decreasing 

(with fixed T"*") increases the viscosity of lubricant. Lubricants with high viscosity has 
higher thresholds for thermal deformation and hence are associated with larger load carrying 
capacity. 

In Figure 3.10 drag force in the sliding direction for different values of k at different inlet 
pressures ranging from = 0 to 8 along the slider have been drawn. From this Figure, 
one can conclude that as k increases the drag force continuously increases for a fixed inlet 
pressure (Pinkus [31]). 



0.4 0.5 0.6 0.7 0.8 0.9 1 

k 


Figure 3.10: Variation of drag force with k at different inlet pressures 

Temperatures along the pad with different boundary conditions have been plotted and 
are shown in Figure 3.11 One can observe that the full flux pad temperatures are higher 
than the corresponding half flux pad temperatures everywhere slightly away from the outlet. 
More prominently, one can notice the existence of thermal boundary layer near the inlet in 
both the cases. In Table 3.3, load carrying capacity for the cases corresponding to no flux 
and half flux models have been presented. Here the no flux model has been analyzed by 
fixing T+ = Tp'*' = 1.5, = I, and = 1 for varying values of k. And the half flux 
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model has been analyzed by fixing = 0.2, T+ = 1.5, = 1 and pt = 1 for different 

values of k. 



Figure 3.11: Pad temperatures at different flux boundary conditions for the case k = 0.4, = 1 


In Table 3.3, load carrying capacity corresponding to above two cases with = 1.5 
and = 0.0 have been presented. Prom these one can notice that both the models 
corresponding to Table 3.3 thermally boost the load carrying capacity of the slider bearing 
with the non zero half flux model having an upper hand. 

Table 3.3: Effect of non zero flux boundary conditions on load carrying capacity at: = 

1.0, Tp+ = T+ = 1.5, pt^l.O ^ ^ ^ 


k 

0.4 

0.6 

0.8 

1.0 

W+(no flux) 

0.7152418 

0.6274425 

0.5569304 

0.500002 

lP+(flux = 0.2) 

0.7301611 

0.6446275 

0.5657586 

0.500002 


The increase in load carrying capacity brought in models corresponding to Table 3.3 is 
due to increase in the viscosity of the lubricant. Prom equation (3.10), one can see that an 
increase in may be possible by either decreasing or . Here as /5‘*‘ is fixed, the 

increase in should be due to decrease in T+. And this is precisely the case as the 
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averaged T of the lubricant from the calculations turn out to be around 1.2. But in the 
simulations corresponding to Table 3.1, the average T+ turns out to be 1.5. 

The study of the effect of considering slider and pad at different temperatures (T/ > T+ 
^ load carrying capacity and frictional drag force of a slider bearing is of 

high practical relevance (ref. Pinkus [31], page 101). Zienkiewicz [42] in his investigation 
considered slider and pad at different temperatures and analytically investigated the issue 
of load carrying capacity with a linear model. As his focus was on parallel slider bearing, 
he has fixed k = 1 and = 0. His studies reveal that if pad temperature (7^" ) is greater 
than the slider temperature (T/), then a suctional effect, which may lead to a drastic fall in 
load carrying capacity and hence a collapse in bearing is possible. In reality, the stationary 
surface (pad) is usually hotter than the moving surface (slider) (ref: Pinkus [31]), so in this 
investigation the following two models have been analyzed: (a) T+ —T^ = 1.5 and = 1.7 
(b) = 1.5 and T/ = 1.7 with pf = 1.0 and 0.4 < A: < 1.0. The load carrying 

capacity corresponding to these two models are presented in Tables 3.4 and 3.5. 

Table 3.4: Comparison of load carrying capacity and drag force for slider and pad setting at; = 
1.5, T+ = 1.7, T+ = 1.5, pf = 1 


k 

0.4 

0.6 

0.8 

1.0 

w+ 

0.7097811 

0.6216089 

0.5546017 

0.500002 

F+ 

0.2908928 

0.3604844 

0.4168325 

0.4692305 


Table 3.5: Comparison of load carrying capacity and drag force for slider and pad setting at: = 

1.5, T+ = 1.5, r+ = 1.7, pt = 1 


k 

■ 0.4 

0.6 

0.8 

1.0 

w+ 

0.7097811 

0.6219599 

0.5541165 

0.500002 

F+ 

0.22216296 

0.2711326 

0.2776720 

0.3321671 


Interestingly, when pf = 1.0, both the models are associated with almost same load car- 
rying capacity for 0.4 k ^ 1.0. However the frictional drag associated with the model(a), 
are higher than those associated with model (b). 
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Another important slider bearing case is, considering inlet lubricant temperature to be 
higher than that of both the pad and slider temperatures. Only two test cases are considered, 
(a) Tt = 1-5, T+ = 1.2, T+ = 1.4 and (b) 7;+ = 1.5, T+ = 1.4, T+ = 1.2. These results 
are presented in Tables 3.6 and 3.7. 

Table 3.6: Vaxiation of load carrying capacity for the setting at: = 1.5, = 1.0 


k 

0.4 

0.6 

0.8 

1.0 

T+ = 1.2, T+ = 1.4 
W+ 

0.7217854 

0.6260851 

0.5550337 

0.5000006 

r+ = 1.4, r+ = 1.4 
1F+ 

0.7198654 

0.6243347 

0.5546634 

0.5000005 


Table 3.7: Influence of non zero flux boundary condition on load carrying capacity for the setting 
at: = 1.5, pf = 1.0 with flux = 0.2 


k 

0.4 

0.6 

0.8 

1.0 

T+ = 1.4 
w+ 

0.7219934 

0.6263183 

0.5552484 

0.5000004 

T+ = 1.2 
PF+ 

0.7165319 

0.6245772 

0.5548232 

0.4999996 


It is observed from Tables 3.6 and 3.7 (with non zero flux) when the slider is at a higher 
temperature as compared to the pad temperature, the load carrying capacity is higher for 
all values of k. Similar result were observed in Tables 3.4 and 3.5. However the cooling of 
the pad and slider (i.e. T+ < and T/ < I)+) yields a higher load as compared to the 
case when T/ > and >T^- 

Other important case is when the inlet lubricant temperature is lower than that of both 
the slider as well as the pad temperature. In this investigation, the following test cases have 
been analyzed (a) — 1.0, = 1-2, T+ = 1.5 (b) = 1.0, T+ = 1.2, — 1.5 

and (c) = 1.5, T/ — = 1-5. The load carrying capacity associated with these 

models for 0.4 < A; < 1-0 and pf = 1-0 have been computed and compared in Figure 3.12. 

The calculated average temperatures of the lubricant with above three models is in the 
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Figure 3.12: Effect of load carrying capacity with k for different slider/pad temperatures at pf = 1 

order: 1.0 < < 1.5 which leads to a reverse order in viscosities 

This justifies the order of the load carrying capacity as seen in Figure 3.12. 

To further understand the effect of zero/nonzero flux boundary conditions on load car- 
rying capacity, simulations have been carried out with pf = 1, T/ = 1.5, % = 1, = 

0, —2.0, 0.2 and the corresponding results are presented in Table 3.8 and Figure 3.13. 

Table 3.8: Effect of non zero flux boundary conditions on load carrying capacity at: = 

1.0, T+ = 1.5,p+ = 1 


k 

0.4 

0.6 

0.8 

non zero flux = 0.2 

0.7301611 

0.6446275 

0.567586 

non zero flux = -0.2 

0.7467269 

0.6854924 

0.6181974 


Prom the plots in Figure 3.13 one may notice that the pad temperatures, at any flux 
namely, positive (influx), zero (adiabatic), and negative (out flux) boundary prescriptions 
are in the decreasing order of magnitude. Also a boundary layer is noticed near the inlet 
boundary. The load carrying capacity in Table 3.8 indicate that the out flux boundary 
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prescription , of the order of - 0.2, with will enhance the load carrying capacity 

in the lubricant and hence constitutes a thermal boosting model. This happens because the 
pad temperatures, as shown in Figure 3.13 are linear and hence the viscosity will be higher. 



Figure 3.13: Effect of flux boundary conditions on pad temperatures at k = 0.4 and pf = I 

For k = 0.4, the streamline pattern corresponding to the inlet pressures pf = 
0, 1, and 5 are presented in Figure 3.14. Prom these patterns one can see that the laminar- 
ity of the lubricant flow is fully lost as the inlet pressure is raised from 0 to 5. The circular 
streamline pattern is observed at pf = 1 can be directly attributed to the bolus like effect 
caused by the higher gage pressure of the order 10^ N/ni^. As the inlet pressure is further 
raised to 5, a bolus like effect with turbulent looking pattern manifests in the lubricant flow 

domain. 

Streamlines corresponding to the cases (a) Tf- = 1.2, T+ = 1.4, T+ = 1.0, k = 0.4 
and pf = 1 (b) Tt = 1.2, Tf = 1.0, Tf = 1.4, k = 0.4 and pf = 1 are shown in 
Figure 3.15. The flow is clearly nonlinear with a bolus like pattern close to either pad or 
slider, on the converging part of the domain, depending on whether T+ > Tf or Tf < Tf. 
These features in (Figure 3.15) the flow pattern can be attributed to the combined effect of 
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(b) 


Figure 3.15: Constant u-velocities for k = 0.4, pf = 1, T+ = 1.2 (a) T+ = 1.4, T+ = 1 (b) T+ = 
1.0, T+ = 1.4 

To get a better understanding of the thermal characteristics of the lubricant, the 
isotherms corresponding to the above two cases are presented in Figure 3.16. Presence 
of the thermal boundary layers on the inlet starting from pad and slider, support the ob- 
servation corresponding to the temperature distribution on the pad. These boundary layers 
are attributed to the difference in the inflow temperature {T^) from that of the pad/slider 
temperatures. 

In Figure 3.17 streamlines and in Figure 3.18 isotherms corresponding to models with 
unequal slider, pad temperatures for fc = 1 (parallel slider bearing) are presented. Again as 
in the tilted pad case one may observe nonlinear and bolus like flow pattern in streamlines 
and boundary layer features in isotherms. One may also observe that the change in slider 
geometry also has its share of influence on the streamline and isotherm pattern. 
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(a) (b) 

Figure 3.18: Toniperaturc contours at A: = 1, = 1, = 1.2 (a) r+ = 1.4, T+ = 1.0. (b) T+ = 

1.0, T+ = 1.4 

In this Chapter, we have chosen some ad-hoc values for the non-zero flux boundary- 
condition on the pad in order to simulate the results (i.e. we have taken 0.2 and also - 0.2). In 
the next Chapter, we have incorporated the heat conduction in the pad. In that simulations, 
the pad flux values turns out to be approximately 0.003. So comparisons are made for load 
carrying capacity with these three values for the setting k = 0.4, pf = 1,T^ = 1,T+ = 1.5 
for Casel and are shown in Figure 3.19. 
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Figure 3.19: Coinpaiison of load carrying capacity with k for various flux prescriptions on the pad 

3.5 Conclusions 

Thermal effects on load carrying capacity in slider bearings has been numerically analyzed by 
SUPGFEM. A thorough investigation of several models has been done. An orderly analysis 
identifying the influence of various factors viz. slider bearing geometry, thermal boost- 
ing, inlet pressure, inlet temperature and temperature dependent viscosity coefficient of the 
lubricant, on the pressure generating mechanism in slider bearings has been provided. Clear- 
ly, the study reveals that the consideration of (a) Lubricant with low values of temperature 
viscosity coefficient (b) Non zero inlet pressure and (c) Tj < < Tp (or) Tj < Tg with 

non-zero flux on pad constitutes a desirable setup for enhancing the Load Carrying Capacity 
of a Slider Bearing. The numerical formulation presented here is suitable for thermal anal- 
ysis of both, high and low speed lubrication flow problems of industrial significance, where 

higher accuracy is required. 




Chapter 4 


Performance of a Tilted Pad Slider 
Bearing Considering Heat Conduction 
in the Pad 


4.1 Introduction 

Lubrication under high load and high speed is common in engineering practice today. In the 
thin film that separates the two rubbing surfaces a large amount of heat may be generated. 
Yet, beyond the classical theory for isothermal and rigid hydrodynamic lubrication, little 
is known about the combined effect of temperature and heat conduction in the solid pad. 
Undoubtly there is a need for a fundamental study of the conditions which surround the 
fluid film (Khonsari and Wang [72]). Several researchers, as discussed in Chapter 3 (page: 
53, 54) have numerically analyzed the thermohydrodynamic performance of slider bearings. 
Studies of Schumack [19] were focussed on the applicability of spectral element scheme for 
thcnnohydrodynamic lubrication (THDL) analysis. Whereas Rodkiewicz and Sinha [58], 
Kumar et al. [73] analyzed the load carrying mechanism of a slider bearings under thermal 
influence by solving nonlinear coupled partial differential equations governing the fluid film 
with either isothermal or adiabatic boundaries or with boundaries having non-zero thermal 
flux prescriptions. 
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Actual bearing configurations have solid pad/slider. To analyze this configuration one 
has to couple the p.d.e s governing the fluid film in the slider bearing with the heat conduction 
equation for the solid pad/slider. Ezzat and Rohde [15] and later Gero and Ettles [18] have 
considered such models for THD lubrication studies. Ezzat and Rohde considered viscosity 
related to two specified temperatures and treated density to be constant. Both viscosity 
and density were treated to be constants by Gero and Ettles [18]. The influence of thermal 
boosting and various geometric configurations on load carrying capacity, drag force etc., of 
the lubricant has been investigated by them. 

In the present Chapter, we investigate the hydrodynamic lubrication of a^nite width 
slider bearings, taking into account not only the temperature variation in the fluid film but 
also the effect of heat conduction to the pad. Influence of thermal boosting, various geometric 
configurations and several practical thermal boundary prescriptions on load carrying capacity 
and frictional drag has been analyzed. Following Brooks and Hughes [68] and Grygiel and 
Tanguy [71] finite element numerical scheme with Streamline Upwind Petrov-Galerkin weight 
functions have been used to preclude the numerical oscillations arising in the flow field. 


4.2 Governing Equations 


The geometry and the coordinate system for a slider bearing has been shown in Figure 3.1. 
Introducing the following non dimensional variables Schumack [19]: 


- 2L 

- 

, X 
■ Ta 

JL 

A+ = XTa, 


+ - 


— sM_ 
^ Uho’ 


,+ ^ 
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(subscript “a” indicates ambient condition). 
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p+ = (p is the gage pressure), 
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kB 
/ + 


X 


V 


IL 

hi 


(where hi is the thickness of the pad). 
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Under the lubrication assumptions mentioned in Section 3.2 and using the above non- 
dimensionalization scheme, the standard Navier-Stokes equations together with the con- 
tinuity equation reduce to the following form: 

_ a , , du^ , 

dx+ dy^ ^ ^ 


^ ' p+ u^) + ^ {p+y+) =Q 

dy+ ^ 


dx+ 


(4.2) 


Combining equations (4.1) and (4.2) with the following boundary conditions 

= 1, v'^ = 0 at = 0 and «+ = 0, v'*' = 0 at j/+ = will lead to the following 
steady state Generalized Reynolds equation for the pressure distribution. 


d h+^ dh+ 

V ..4- u ) 


dx'^ 6 (9x+ dx'^ 

with the boundary conditions 

—pf at = 0 and p'*' = 0 at a:'*' = 1 


(4.3) 


(4.4) 


The energy equation is as follows 

+ , +3T+ , ^dT+ , 1 d^T+ , PrE, + , eu+ ,, , P, E. dp+ ,, ,, 

(4.5) 

where 


p+ = 1 - A+ ( T+ - 1 ), = exp [ -,0+ ( T+ - 1 ) 


(4.6) 


Equation governing the heat conduction in the pad is 




dx 


+2 


dy' 


+2 


(4.7) 


In view of the various industrial applications, the energy equation has been solved under 
the following different temperature boundary conditions (see Figure 3.1 for AB, BC, AD, 
CD, DF, CE, EF positions): 
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• Casel: 

Lubricant boundary conditions: 

T+ = T+ = 1.3 ( or 1.2 or 1.1 or 1.0 ), on AB 
1 t = 0, onBC 
T+ = I^+ = 1.3, on AD 
Pad boundary conditions: 

= 1.3,, on DF 
T+ = 1.3„ onCE 
= 1.3„ on EF 

Interface boundary conditions: 

T+ = V and ks^ = kf^, onCD 

• Case2: 

Lubricant boundary conditions: 

T+ = T+ = 1.3 ( or 1.2 or 1.1 or 1.0 ), on AB 
gi = 0, on BC 
T+ = = 1.3, on AD 

Pad boundary conditions: 

^ = 0, onDF 
1^ = 0, CE 
Tp+ = 1.3, on EF 
Interface boundary conditions: 

T+=V “d = onCD 

• Case3: 

Lubricant boundary conditions: 


T+ = T+ = 1.3 ( or 1.2 or 1.1 or 1.0 ), on AB 


dT+ 


= 0 , 


on BC 


T+ = Tt = l.3, on AD 
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Pad boundary conditions: 
?^ = 0, onDF 


^ = 0, onCE 
¥ = 0, onEF 


Interface boundary conditions: 

T+ = Tp+ and k,^ = kf^, on CD 


• Case4: 


Lubricant boundary conditions: 

T+ = T+ = 1.3 ( or 1.2 or 1.1 or 1.0 ), on AB 
gi = 0, onBC 
T+ =Tt = 1.3, on AD 
Pad boundary conditions: 

^ = 0, onDF 
gj = 0, onCE 

-ks^ = -H (r/-l), onEF 


Interface boundary conditions: 


Tf = T* and K ^ = k, ^ 




The load carrying capacity (PF"*") and the frictional drag have been determined using 
the equations (3.11) and (3.12). 


4.3 Finite Element Formulation and Solution Proce- 
dure 


Let 0 the domain of interest be bounded in with a piecewise smooth boundary F. Let 
(a;+, y+) denote the vector of spatial coordinate of a generic point in Q. and n is the outward 
unit normal vector drawn at any point (a:+, y"*") on F. Now the variational principle 11 (p) for 
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the Reynolds equation (4.3) in is given by; 


n(p) = f I 

JQ 




dp 


,+ 




(4.8) 


The weighted residual formulation of the energy equation is given by: 


/„ 


n ox^ 


+ v 


dT+ 1 d'^T+ PrE.. 
dy+ ^ Pe ay+2 


+ r ^ dW^ ^2 Pr E, 


U+ ( — ) 


r 4-/C 4- 


U ' 


dx'^ 


dn = o 

(4.9) 

where, is the modified SUPG weight function. It is to be noted that are 

dependent on T+ and are defined in equation (4.6). 

And the Galerkin formulation for the heat conduction in the pad is given by: 

= o (4.:0) 

Now to reduce the lubrication problem to one of finite number of unknowns, the solution 
domain is divided into bilinear rectangular elements and the field variables (pressures, ve- 
locities and temperatures) are expressed in terms of approximating trial functions Ni within 
each element. Now Q is discretized into Ne bilinear rectangular elements such that: 

U^i a^ = Q O® = 0 (4.11) 


where, denotes the interior domain of an element. Let F® be the boundary of f2®. The 


discretized representation of the field variables are given by: 

N,i 




E 


u 


k ’ 


V* Si 


E Pi 


A:=l 


A:=l 


Nel 

E Jv. pt\ 

k=l 


Nel 

r+ Si Y, Nc T^'. 

k=l 


Tt « 


E 

k=l 


where, N^i = 2 for pressure case and Ngi = 4 for velocities and temperatures. Introducing 
these elemental discretization details into equations (4. 8)- (4. 10) we obtain element level 


representations. The element level representation of equation (4.8) is given by: 

I < iS 1 ^ ( I ^ ( E .r ) } ^«=o (4.12) 
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on minimizing the equation with respect to variables p® and on further simplification one 

dNf 


would arrive at the following equation. 

/• dm dNf 


^ 6p+ dx+ dx+ ine dx+ 


(4.13) 


Now with the SUPG analysis of Brooks and Hughes [68] and Grygiel and Tanguy [71] one 
would obtain the following discretized elemental energy equation. 


J (EKTf - 1 )) [ ^ ^ 

j k k k 


+ 


1 , dwi dNf,, f ,-PrE, 


Pe dy+ dy+ 


)} prdn^ = JjW-^iexT>{-p-^(£NfT,^‘ - l)))(^(E^^f )) 


dx+ dy+ 


d 


- CZKO in' = 0 

k k 


(4.14) 


where, 


ic 'll W-- 

W = Wi + Pi Wi = Ni and Pi — — ;" ■ ■ ■ . — 

u 


And the upwind parameter k is calculated using elemental dimensions and elemental velocity 
(Grygiel and Tanguy [71] and Brooks and Hughes [68]. 

And the discretized elemental level heat conduction equation is given by: 

(4.15) 


i 


r ,,+.2 dwi dNj‘ dwi aiV/ we 


Now based on equations (4.13) - (4.15) numerical simulations have been carried out on 
a 30 X 30 finite element mesh system. The results have been obtained to an accuracy of 
e = 5.0 X lO-'^ i.e 

I ^New _ ^Old I 

{ ' )) 

\ fj \ 

for j = 1, 2, • • • , n, ■ ■ where 7j represents any of the field variables p+, 
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4.4 Results and Discussion 

The simulations have been carried out on 30 x 30 mesh system. Computer generated plots 
are obtained for various values of k, pf and with the fixed values for the parameters Pr 
(Prandtl number) and Ec (Eckert number). The results have been analyzed for load carrying 
capacity, drag force and velocity, pressure, temperature fields and are presented in the form 
of plots and Tables. 

The non dimensional pressures for Casel at different locations in the sliding direction 
with k = 0.4, = 0 for different slider temperatures (T^"^) are shown in Figure 4.1. 



Figure 4.1: Pressure distribution along the sliding direction of the bearing at pd = 1, ft = 0.4 

From this Figure one can notice that for fixed values of a;'*', lowering the slider tempera- 
ture {Ts^) below the inlet lubricant temperature (T)+) enhances the pressure. So maintaining 
the slider at temperature lower than that of inlet lubricant temperature leads to thermal 
boosting of the slider bearing performance. In Figure 4.1 one can notice that when the non 
dimensional slider temperature is kept at = 1.0 and the inlet lubricant temperature at 
= 1.3, peak pressures are observed at a;+ = 0.8. Both of these observations can be 
attributed to the fall in the average temperature of the lubricant in the domain under con- 
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sideration, on lowering the slider temperature. Prom the expression for viscosity (;Li+(T+)) 
it is clear that for fixed /?+ lowering T+ leads to higher values for viscosity. Higher the 
viscosity higher will be the threshold for thermal deformation of the lubricant and this leads 
to improved performance. Such a fail in the lubricant temperature can be ascertained by 
tracing the temperature distribution in the form of isotherms of the lubricant for different 
slider temperature set ups. So in Figure 4.2 isotherms for the set ups T+ = 1, 1.1, 1.2, 1.3 
with T+ = 1.3 for all these cases are presented. 



(c) (d) 

Figure 4.2: Isotherms at (a) = 1.0, (b) r+ = 1.1, (c) T+ = 1.2, (d) = 1.3 
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These plots fully support the observations regarding the lubricant temperature fall. 
One can notice that around x'^ = 0.8, 0 < y'^ < for all slider temperature set ups, 
average temperature distribution reaches a minimum and this lower bound can be lowered 
by decreasing T,'*' which in turn can aid the pressure build up. Also on lowering Ts^ a 
thermal boundary layer is seen to develop adjacent to the inlet. In the region R between 

0. 2 < a; < 0.5, 0.35 < y < 0.5 in the Figure 4.2 circular isotherms with increasing inner 
temperatures are noticed. In one case Figure 4.2 (d) the temperature goes beyond inlet 
temperature 1.3. These are hot spots in the domain and they indicate the presence of small 
recirculating fluid zones. To verify this, u-velocity contours ( iso u-velocities or streamlines 
as V « 1 ) are presented in Figures 4.3 and 4.4 for = 1.3, T+ = 1.0, 1.1, 1.2, 1.3, pf = 

1, A: = 0.4, = 0. 




(a) 


(b) 


Figure 4.3: u-velocity contours at (a) T+ = 1.0, (b) T/" = 1.1 
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Figure 4.4: u-velocity contours at (c) T/ = 1.2, (d) = 1.3 

In these Figures one can notice the naanifestation or positive indication of recirculation 
zones in the region R. These circulation zones get prominent as TJ' is set closer to The 
isotherms in Figure 4.2 indicate that as the recirculation zones in region R in Figures 4.3 and 
4.4 get prominent, the hot spots or circular isotherms too get stronger. From the Figure 4.2 
(a) and (b) corresponding to = 1.0 and 1.1 one can notice that in the inlet region where a 
thermal boundary layer with plume like isotherm manifests, flow field gets concentrated with 
fluid boundary layer immediately followed by a long and slender recirculation zone. This 
can be attributed primarily to the thermal boundary conditions and partly to the geometric 
configuration of the slider bearing. In Figures 4.3 and 4.4 one can also see that as is 
raised boundary layer adjacent to the inlet vanishes and flow in the domain gets accelerated. 
The flow pattern for various- values of T+ suggest that on lowering T/ flow gets more laminar 
excepting close to the inlet. 

In Figure 4.5 the frictional drag force (F’+) corresponding to the above simulations 
along the sliding direction have been plotted. It is clear that T/ = 1.0, = 1.3 not only 

leads to higher pressure generation but also to reduced frictional drag. From the formula for 
frictional drag force given in equation (3.12) as /j,'^ is already seen to increase on lowering 
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slider temperature any reduction in F'^ should be due to fall in the velocity gradient . 
The u-velocity contours in Figures 4.3 and 4.4 show that such a fall in velocity gradient does 
indeed take place. 



Figure 4.5: Frictional drag force along the sliding direction of the bearing at pf = = 1-3, k — 

0.4, Casel 

The combined effect of thermal boosting and geometric configuration of slider bearing on 
load carrying capacity and the frictional drag are presented in Figures 4.6 and 4.7 respectively 
for Case4 at = 1. These plots show that decreasing either k or T/ or both k and T+ 
simultaneously enhances the load carrying capacity but lowers the frictional drag values. 
Further the combined influence of k and T+ on frictional drag and load carrying capacity is 
stronger than their segregated effects. While the combined influence of k and on frictional 
drag is clearly seen at all k, it is predominantly seen for OA < k < 0.6 on load carrying 
capacity. 
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Figure 4.7: Variation of frictional drag force with k at = 1, = 1.3, Case4 
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In Figure 4.8 the combined influence of thermal boosting and the geometry of the slider 
bearing on load carrying capacity for pf = 0 is presented. One may notice that these results 
are similar to those obtained at = 1 (see Figure 4.6). Further these are in conformity 
with the remarks made by Pinkus and Sternlicht [74]. 



Figure 4.8: Pressures along the pad at = 0, fc = 0.4 


Also from Figures 4.6 and 4.8 it is seen that the thermal boosting has greater influence 
on load carrying capacity for smaller inlet pressures. In the case k = 0.4, thermal boosting 
i.e = 1.1, has brought in nearly four fold increase in the load carrying capacity. It may 
also be interesting to trace the pressures along the sliding direction on the pad for k = 0.4 
at pf = 0. 

From Figure 4.8 one may notice that for all T+ prescriptions maximum pressure on 
the pad manifest around x'^ = 0.75 to 0.8 which is close to the outlet. This result is very 
much in conformity with the standard results reported on load carrying capacity for slider 
bearings (Rodkiewicz and Sinha [58]). A more interesting feature which one would notice in 
Figure 4.8 is that the thermal boosting effect on the pressure along the pad is also maximum 
around x'^ = 0.75 to 0.8. 

Load carrying capacity of parallel slider bearing (k = 1) with zero inlet pressure {pf ) has 
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lot of significance in several applications. So the influence of thermal boosting on pressure 
generation in the fluid film for A: = 1 with pf = 1 and Case4 boundary conditions has been 
analyzed and presented in Figure 4.9 for various values of T+. 



Figure 4.9: Load carrying capacity at different slider temperatures with pf = 0,k = l, Cased 

It is clear that on reducing T+ from 1.4 to 1.0, the load carrying capacity is boosted 
more than two and half times, i.e from 3 x 6“°® to 7.75 x 6"°® . This result supports the 
view that pressure generation is very much possible in parallel slider bearings. Here it is to 
be noticed that in this study, p"*" are functions of temperature. 

Next the influence of solid pads with various thermal boundary conditions as described 
by Casel ( pad with isothermal side and isothermal top face ), Case2 ( pad with adiabatic 
sides and top ) and Cased ( pad with adiabatic sides and exposed top ) on the load carrying 
capacity and frictional drag force in the fluid film has been analyzed for = 1.3, 1.2, Tj' — 
1 . 3 , p+ = 1, k = 0.4, 1. These results are presented in Tables 4.1 - 4.4. 
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Table 4.1: Comparison of load carrying capacity and frictional drag force at different boundary 
conditions for ~ 1.3, r+ = l.Z,pf = 1,A: = 0.4 


Boundary conditions 

Casel 

Case2 

Case3 

Case4 

w+ 

0.7364519 

0.7364491 

0.7364560 

0.7364571 

F+ 


0.2086195 

0.2086191 

0.2086179 


Table 4.2: Comparison of load carrying capacity and frictional drag force at different boundary 
conditions for = 1.3, T+ = 1-3, = 1,A: = 1 


Boundaxy conditions 

BIIIIIIISSSSH 



Case4 


lll!^ 



0.5000002 , 

F+ 

0.4217472 

0.4217471 

0.4217361 

0.4217463 


Table 4.3: Comparison of load carrying capacity and frictional drag force at different boimdary 
conditions for = 1.3, T+ = 1.2, pt = 1,A; = 0.4 


Boundary conditions 

Casel 

Case2 

Case3 

Case4 

w+ 

0.7599708 

0.7599646 

0.7588058 

0.7587989 

F+ 

0.866119 

0.0865929 

0.08708 

0.08708 


Table 4.4: Comparison of load carrying capacity and frictional drag force at different boundary 
conditions for = 1.3, T/ = 1.2,p^ = l,fc = 1 


Boundary conditions 

Casel 

Case2 

Case3 

Case4 

w+ 

0.4999992 

0.5000011 

0.5000002 

0.5000002 

F+ 

0.3655055 

0.3654333 

0.3654611 

0.3654560 
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Prom Tables 4.1 - 4.4, both at — 1.3 and = 1.2 only a marginal variation in load 
carrying capacity and frictional drag force is noticed. To further investigate this average 
temperature flux across the solid pad and lubricant interface are evaluated for Casel , 
Case2 , Case3 and Case4 setups and these are found to be in the range of 0.002 to 0.005. To 
assess the influence of flux boundary values naturally getting set on the interface, a special 
slider bearing conflguration without solid pad has been considered and a temperature flux 
boundary condition is prescribed on the pad side. 

Experiments are carried out for different temperature flux prescriptions setting the 
remaining parameters at: pt = 1, = 1-3, T+ = 1.2 and 0.4 < A; < 1. Results in the 

form of load carrying capacity obtained for flux prescriptions of 0.002 to 0.004 are presented 
in Table 4.5. One can see that small variation in flux values from 0.002 to 0.004 does not 
bring in a visible change in the load carrying capacity. These results agree with Rohde and 
Kong [49]. 


Table 4.5: Comparison of load carrying capacity at different flux values and at different values of 
k with pf = 1. 


Flux 

k = 0.4, W+ 

k = 0.6, W+ 

A: = 0.8,1F+ 

k = l,W+ 

0.002 

0.778611 

0.6421080 

0.5604001 

0.500002 

0.003 

0.778603 

0.6421075 

0.5603997 

0.500002 

0.004 

0.778595 

0.6421067 

0.560389 

0.49999 


Isotherms corresponding to the slider bearing with solid pad for the Casel and Cased 
are presented in Figure 4.10 (a) and (b). The isotherm patterns close to the interface clearly 
show the temperatures in the vicinity of interface, both on pad and lubricant side have only 
marginal changes and hence give raise to small flux values. Further it must be noticed that 
on the pad side little above the interface multiple circular isotherm pattern with almost 
fixed isotherm values are present. These go to say that fluxes on the interface, excepting 
close to the inlet and outlet will be very small. However, the influence of the Cased on 
the temperature distribution in pad may be noticed from the isotherms density close to the 
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adiabatic sides and also from their magnitudes in the interior of the pad. 




Figure 4.10: Isotherms at = 1.3, A: = 0.4, = 1 for (a) Casel, T+ = 1.3 and (b) Case4, T+ = 
1 . 2 . 

In Table 4.6, the effective load caxrying capacity for the case with conduction and no 
conduction in the pad are shown by fixing — 1,TJ'+ — 1-5, pf = 1, Case4, flux = 0.003. 
From this Table 4.6, we notice that the load carrying capacity for the case of no conduction 
in the pad are higher for values of k = 0.4, 0.6, 0.8. Whereas for parallel case {k = 1) no 
change in the load carrying capacity is observed. 


Table 4.6: Comparison of load carrying capacity for the case with conduction and no conduction 
in the pad for pf — l-iTf = l,Tf = 1.5, flux = 0.003, Case4 



0 

II 

k = 0.6 

00 

0 

II 

k = 1 

W'^ (with conduction) 

0.7365377 

0.6255385 

0.5560876 

0.500002 

(no conduction) 

0.7786063 

0.6421075 

0.5603997 

0.500002 
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4.5 Conclusions 

Load carrying capacity and frictional drag of a high speed slider bearing with a solid pad 
for various thermal boundary conditions in conjunction with thermal boosting for various 
geometric configurations have been numerically analyzed using FEM. Density and the Vis- 
cosity in the mathematical model describing the slider bearing configuration are treated as 
functions of temperature. SUPG weight functions are found to be effective in precluding 
the numerical oscillations and in obtaining a convergent solution. Load carrying capacity of 
slider bearing can be thermally boosted by setting slider at temperatures lower than that of 
incoming lubricant. Thermal boundary layer in the inlet region, hot spots and recirculation 
zones in region R ( 0.2 < x'^ < 0.5, 0.35 < y"*" < 0.5 ) sensitive to are noticed in the 
temperature and fiow fields respectively. Lowering of slider temperature takes the fiow closer 
to the laminar regime and can make slider bearings more durable. Though solid pads with 
different boundary settings as chosen in this study from application point of view alters the 
the temperature field in the pad the thermal flux across the fluid-solid interface changes 
marginally and thus the load carrying capacity of the slider bearing remains unaffected.. It 
is to be noted that the isothermal pads employed in the present study have been set at 
temperature similar to that of inlet flow. 



Chapter 5 

Thermohydrodynamic Analysis of a 
Tilted Pad Slider Bearing Considering 
Heat Conduction in the Pad and 
Slider^ 


5.1 Introduction 

Tilted pads are widely used in many types of machinery because of their excellent stabili- 
ty, superior durability and high load carrying capacity. Numerous work has been done to 
investigate the behavior of tilted pads from various points of view. It is believed that perfor- 
mance of a tilted pad is affected by many factors, such as pad conduction, heat conduction 
in the moving slider, flux boundary conditions, effect of the inlet pressure buildup, etc. 

The importance of thermohydrodynamic lubrication of slider bearings has been high- 
lighted in the previous Chapter. In most of the works cited earlier the heat conduction 
in the slider and pad has not been considered. Dowson and Hudson [44] re-examined the 
mechanism of slider bearings allowing for the heat transfer to the bearing solids to deter- 
mine the fluid film thermal boundaries. More recently Ezzat and Rohde [15] analyzed the 
thermohydro dynamic behavior of a steadily loaded inclined plane slider, taking into account 


^Accepted for presentation and publication in the STLE Annual Meeting Proceedings, May 2000 
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the finite width of the bearing and heat transfer to the solids. Kumar et al. [73] have nu- 
merically analyzed the thermohydrodynamic performance of a slider bearing under thermal 
influence by solving nonlinear coupled partial differential equations governing the fluid film 
with either isothermal or adiabatic boundaries or with boundaries having non zero thermal 
flux prescriptions. Later Kumar et al [73] also analyzed the affect of heat conduction in 
the pad on bearing characteristics. But actual bearing configurations have slider and pad 
which are solids of finite thickness.To analyze this configuration one has to couple the p.d.e’s 
governing the fluid film in the slider bearing with the heat conduction in the pad and as well 
as in the moving slider. 

Rohde and Kong [49] and later Gero and Ettles [18] have considered only heat conduc- 
tion in the pad but not both to the pad and slider. Rohde and Kong [49] have assumed 
viscosity related to two specified temperatures and density to be a constant, whereas Gero 
and Ettles [18] considered viscosity and density to be constants. The influence of thermal 
boosting and various inclinations of pad configuration on load carrying capacity, drag force 

etc., of the lubricant has not been investigated by them. 

In the present Chapter, temperature sensitive density and viscosity, which brings in 

extra pseudo compressibility terms into the energy equation have been analyzed. Influence 
of thermal boosting, various geometric configurations and several realistic thermal boundary 
prescriptions on load carrying capacity and frictional drag has been analyzed. Following 
Brooks and Hughes [68] and Grygiel and Tanguy [71] finite element numerical scheme with 
Streamline Upwind Petrov-Galerkin weight functions has been used to preclude the numerical 
oscillations and results are depicted through computer generated plots and tables. 

5.2 Problem Formulation 

The geometry and the coordinate system for a slider bearing analyzed in this Chapter is 
shown in the Figure 5.1. 

Introducing the following non dimensional variables Schumack [19]: 



5.2 Problem Formulation 


100 



Figure 5.1: Geometry of the slider bearing considering heat conduction in the pad and moving 
slider 
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hf = ^. (where hi is the thickness of the pad). 
hfi — (where hn is the thickness of the slider). 

Under the assumptions given in Chapter 3 (section 3.2) and using the above non- 
dimensional scheme, the standard Naiver-Stokes equations together with the continuity e- 


quation reduce to the following form: 

dp'^ 


9x+ dy+ 


d , , du^ . 


(5.1) 




(5.2) 
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Combining equations (5.1) and (5.2) with the following boundary conditions 
«+ = 1, = 0 at ?/+ = 0 and = 0, = 0 at 

will lead to the following steady state Generalized Reynolds equation for the pressure distri- 
bution. 


d . dp^ , _ dh^ 
dx'^ 6 dx'^ 


(5.3) 


with the boundary conditions 

=Pt at a:"*" = 0 and p"*" = 0 at a:'*' = 1 (5.4) 


The energy equation is as follows 

+ .+ dT+ + ar+ 1 d^T+ ^PrE, + _ 2 

^ ^ ^ dx+ ^ dy+ ^ Pe dy+^ ^ Pe ^ ^ dy+ ^ ^ Pe 


dy+ 


du^ .9 


r 4 - 




dx+ 


where, and /x''" are given by: 


)+ = 1 - A+ ( T+ - 1 ), = exp [ -,5+ ( T+ - 1 ) ] 


(5.5) 


(5.6) 


Equation governing the heat conduction in the pad is 

, .^,2 (fix* _ 

And the equation governing the heat conduction in the slider is 


., + .2 d^T+ , d‘^T + 

ll) g ^„+2 ^^,,+2 


0 


(5.7) 


(5.8) 


In view of the various industrial applications, the energy equation has been solved under 

the following different temperature boundary conditions: 

• Casel: 

Lubricant boundary conditions: 
fi = 0, onBC 


T+ = r+ = 1.3, on AD 
Pad boundary conditions: 
Tp+ = 1.3, onDF 



5.2 Problem Formulation 


102 


T+ = 1.3, onCE 
2 ;+ = 1.3, onEF 
Slider boundary conditions: 

T+ = 1.3, on AG 

r+ = 1.3, onBH 

T+ = 1.3, onGH 

Interface boundary conditions: 

r+ = r/ and on CD 

T+ = T/ = 1.3 ( or 1.2, or 1.1, or 1.0 ) and k' s ^ = kf on AB 
• Case2: 

Lubricant boundary conditions: 

= 0, on BC 


r+ = i;+ = 1.3, on AD 
Pad boundary conditions: 


= 0, on DF 


dTp+ 
dx 

^= 0 , CE 

r,+ = 1.3. onEF 

Slider boundary conditions: 


= 0, on AG 


err^-^ _ 

dx 

^ - 0, BH 


dx 


T¥ 


T+ = l.3, onGH 
Interface boundary conditions: 

T+ = T+ and k,^ = kf^, on CD 

T+ = T/ = 1.3 ( or 1.2, or 1.1, or 1.0 ) and k' s on AB 

• Case3: 

Lubricant boundary conditions: 
fi = 0, onBC 
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T+ = i;+ = 1.3, on AD 


Pad boundary conditions: 
1^ = 0, onDF 
^ = 0, onCE 
^ = 0, onEF 
Slider boundary conditions: 
g^ = 0, on AG 
= onBH 


®^ = 0, onGH 


Interface boundary conditions: 
r+ = T/ and k,^ = kf^, on CD 

r+ = Ts'^ = 1.3 ( or 1.2 or 1.1 or 1.0 ) and k's = kf “ 


• Cased: 

Lubricant boundary conditions: 
gi = 0, onBC 
r+ = 7^+ = 1.3, on AD 
Pad boundary conditions: 

^ = 0, onDF 

dx^ ’ 

^^^ = 0, onCE 


dy” 


V dy+' 


dx 

-k ^ 

dfi 

Slider boundary conditions: 


H ( V - 1 ), on EF 


dTs+ _ 
dx 




= 0, on AG 


f^ = 0, onBH 

-k’,^ = -H (T/-1), onGH 
Interface boundary conditions: 
r+ = r/ and A:, = on CD 

r+ = Ts'^ = 1-3 ( or 1.2 or 1.1 or 1.0 ) and k' 


dT,+ 

s dy + 


'^f dy+^ 


on AB 


on AB 
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The load carrying capacity {W'^) and the frictional drag (F+) have been determined from 
equations (3.11) and (3.12) respectively. 

5.3 Numerical Method 


Let 0 the domain of interest be bounded in with a piecewise smooth boundary F. Let 
y+) denote the vector of spatial coordinate of a generic point in 0 and n is the outward 
unit normal vector drawn at any point (x'^, y+) on F. Now the variational principle n(p) for 
the Reynolds equation (5.3) in Q is given by: 


12 dx'^ 




dp'^ 

dx-^ 


dO, 


(5.9) 


The weighted residual formulation of the energy equation is given by (Reddi [70]). 

fr7r^+,4.dT+ ^dT+, 1 d^T+ PrBc PrEc 

(^^r) ® (5.10) 






dy+’ Pe dy+^ Pe ^ ^dy+ 

where, Wi is the modified SUPG weight function. It is to be noted that p"*", are 
dependent on T'*' and are defined in equation (5.6). Galerkin formulation for the heat 
conduction in the pad is given by: 

L ^ ay‘*‘ ° 


(5.11) 


And the Galerkin formulation for the heat conduction in the moving slider is 

r , ,^+,2 d^T+ _ 8^%-^ , 

L" ^ ^r"+" dy"^^ ° 


(5.12). 


Now to reduce the lubrication problem to one of finite number of unknowns, the so- 
lution domain is divided into bilinear rectangular elements and the field variables (pres- 
sures,velocities and temperatures) are expressed in terms of approximating trial functions 
Ni within each element. Now 0, fi', are discretized into iVe bilinear rectangular elements 
such that: 


nfii n‘ = 4' 



5.3 Numerical Method 


105 


= Q' nfji fy® = (!> 

sf'" = n" = (^ 


where, 03 denotes the interior domain of an element. Let F® be the boundary of 03. The 


discretized representation of the field variables are given by: 

^cl 


u: 


E “f. 


v'^ ~ 


E N, vf, 


Pt 


k=l 


k=l 


Nel 

E Pf 

k=l 


N'cl 


N,i 


Tt «= E r+‘. « E . 

fc=l A:=l 

where, Na = 2 for pressure case and Ne,i — 4 for velocities and temperatures. Introducing 
these elemental discretization details into equations (5.9) - (5.12), we obtain element level 


representations. The element level representation of equation (5.11) is given by: 

1 ^ ( I "7 ) ^D}^ = 0 


(5.13) 


on minimizing the equation with respect to variables p® and on further simplification one 
would arrive at the following equation. 


V f — ^ ^ 

^ iu' 6 (9a:+ 


pfdO^= [ h+ 


aiv/ 

5rr+ 


dO® 


(5.14) 


Now with the SUPG analysis of Brooks and Higher [68] and Grygiel and Tanguy [71] one 
would obtain the following discretized elemental energy equation. 

E /„.{ W' (1 - A" - D) I (Eufjva ^ + (T^vt-Ni) 1 


+ A' = i ^(<^p(-^*(Eww‘ - 

d 


k k 


(5.15) 
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where, 


Wi = wi + P, Wi = Ni &nd Pi = 




And the upwind parameter k is calculated using elemental dimensions and elemental velocity 
(Grygiel and Tanguy [71] and Brooks and Hughes [68]). The discretized elemental level heat 
conduction equation in the pad is given by: 

f r /7 4-^2 dWi dN-i^ dWi dNP,^.e ,^'e „ / s 

i- ' ^ ^ ^ 

And the discretized elemental level heat conduction equation in the moving slider is 


[». W.) 


2 dWi dNj^ dWi diVj® W'e 

dx"-^ dx''+ dy"^ dy"-^^ 


0 


(5.17) 


Now based on equations (5.14) - (5.17) the numerical simulations have been carried out 
on a 30 X 30 finite element mesh system. The results have been obtained to an accuracy of 
e = 5.0 X 10"^ i.e 

I yyNew __ r^Old I 

(^a^( ' )) < ^ 


I 7f 


for i = 1, 2, • • • , n, • • •. where, 7 i represents any of the field variables Tf. 


5.4 Results and Discussion 


The simulations have been carried out on 30 x 30 mesh system. The results have been 
computed for load carrying capacity, drag force, velocity, pressure and temperature fields. 
These are presented in the form of computer generated plots and tables on a 30 x 30 mesh 
system. 

The non dimensional load carrying capacity and frictional drag force obtained for k = 
0.4,pt = l, 2;+=r+ = i.3 are shown in Table 5.1 for various thermal boundary conditions 
viz. Casel, Case2, Case3, Cased. 

Considering Cased, as a natural heat conduction to the environment, the load carrying 
capacity at different values of and for two different inlet temperatures are shown in 
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Figure 5.2. Prom this Figure one can notice that for fixed values of x'^, lowering the slider 
temperature (T+) below the inlet lubricant temperature lead to thermal boosting of the 
slider bearing performance. 


Table 5.1: Comparison of load carrying capacity and frictional drag force at different thermal 
boundary conditions for pf = 1, k = 0.4, = 1.3 


Boundary conditions 

Casel 

Case2 

CaseS 

Case4 

w+ 

0.7371731 

0.7371749 

0.7371752 

0.7371781 

F+ 

0.2016197 

0.2016298 

0.2011301 

0.2011315 



Figure 5.2: Pressure distribution along the sliding direction of the bearing at A: = 0.4, pf — 1, 
Case4 

In Figure 5.3 pad distribution is shown for different slider temperatures. One can notice 
that when the non dimensional slider temperature is kept at r+ = 1.0 and the inlet lubricant 
temperature at Tf' = 1.3, peak pressures are observed at x'^ = 0.8 for pf = 0. Both these 
observations can be attributed to the fall in the average temperature of the lubricant in the 
domain under consideration, on lowering the slider temperature. Prom the expression for 
viscosity p (T+), it is clear that for fixed lowering T+ leads to higher values for viscosity. 
Higher the viscosity higher will be the threshold for thermal deformation of the lubricant 
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and this leads to improved performance of the slider bearing. 



Figure 5.3: Pressure distribution along the sliding direction of the bearing at A: = 0.4, pf = 0, 
Case4 


The combined effect of thermal boosting and geometric configuration of slider bearing 
on load carrying capacity and frictional drag are presented in Figure 5.4 and 5.5 respectively 
for Cased pf = 1. 

From the plots given in Figure 5.4 and 5.5, it is clear that decreasing either k or T+ 
or both simultaneously enhances the load carrying capacity but lowers the frictional drag 
values. Further the combined influence of k and T+ on frictional drag force is clearly seen 
at all values of k, it is predominantly seen for OA < k < 0.8 on load carrying capacity. In 
the case k = 0.4, thermal boosting i.e T+ = 1.1, has brought in nearly two fold increase in 
the load carrying capacity. It may also be interesting to trace the pressures along the sliding 
direction on the pad for k = 0.4 at pf = 0. Hence in Figure 5.6 the corresponding pressures 
along the pad are presented. 

From the Figure 5.6 one may notice that a reduction in the = 1.2 does not shift 
the location of the maximum pressure on the pad significantly. This result is very much in 
conformity with the standard results for slider bearings (Pinkus and Strenlicht [74]). 
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Figure 5.6: Pad pressures along the sliding direction of the bearing at A: = 0.4, pf =■! 

An interesting feature which one may notice in Figure 5.6 is that in the thermal boosting 
case the pressure peak is more pronounced (around = 0.75 to 0.8). 

The influence of solid pad and moving slider with various thermal boundary conditions 
described by Casel ( pad with isothermal side and top face, slider isothermal side and bottom 
face ), Case2 ( Pad with adiabatic sides and top, slider with adiabatic sides and bottom) and 
Cased ( pad with adiabatic sides and exposed top, slider with adiabatic sides and exposed 
bottom ) on the load carrying capacity and frictional drag force in the fluid film has been 
analyzed for = 1.3, = 1.3, pf = 1, k = 0.4, 1. These results are presented in Tables 

5.2 - 5.4. Table 5.2 depicts the results for parallel slider bearing. It is seen that a finite load is 
generated, due to the assumption of the fore=region pressure. This load remains unaffected 
by the various thermal boundary conditions. It is observed that even if we reduce the slider 
temperature below the inlet temperature there is no change in the load carrying capacity 
(See Tables 5.3 and 5.4), however a slight variation is observed in the frictional drag force. 
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It is observed that even if we reduce the slider temperature below the inlet temperature 
there is no change in the load carrying capacity (See Tables 5.3 and 5.4), however a slight 
variation is observed in the frictional drag force. 


Table 5.2: Comparison of load carrying capacity and frictional drag force at different thermal 
boundary conditions for = 1, A: = 1, = r+ = 1.3 


Boundary conditions 

Casel 

Case2 

Case3 

Case4 

w+ 

0.4999992 

0.4999993 

0.5000002 

0.4999992 

F+ 

0.417500 

0.4176240 

0.40002 

03951091 


Table 5.3: Comparison of load carrying capacity and frictional drag force at different thermal 
boundary conditions with = 1, k = 0.4, = 1.3 for Casel 



k = 0.4 

k = 0.6 

k = 0.8 

k = 1.0 

w+ 

0.785592 

0.6554149 

0.5655332 

0.4999992 


ii.0S2?'I9 

■■().i<)()3;()2 , 

(1.27-.' 12: 

1 T.'ii'i) 

1 1 


Table 5.4: Comparison of load carrying capacity. and frictional drag force at different thermal 
boundary conditions withp^ = 1, A: = 0.4, = 1.2, T/ = 1.1 for Casel 



k = 0.4 

k = 0.6 

k = 0.8 

k = 1.0 

w+ 

0.7855930 

0.6554189 

0.5655632 

0.4999992 

F+ 

0.0862 

0.19312 

0.28127 

0.44700 


To further investigate this, average temperature flux across the slider-lubricant and 
pad-lubricant are evaluated for the Case4 thermal boundary condition. These interface 
temperatures at various sliding directions are shown in the Figure 5.7 and 5.8. 
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Figure 5.8; Pad-lubricant interface temperature gradients along the sliding direction of the bearing 
a,t k = 0.4, pf = l 
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From Figure 5.7 one can notice that there is a very little variation in slider-lubricant 
interface temperature gradients, thus the thermal boosting is not sufficient to bring forth 
any significant changes in the load carrying capacity. Whereas in Figure 5.8 absolutely no 
changes in the pad-lubricant interface temperature gradients are observed. 

To further assess the influence of flux boundary conditions, pad temperatures at various 
values of x+ are shown in Figure 5.9. Prom this Figure one can conclude that although 
the pad temperatures are slightly high for = 1.3, T+ = 1.3 compared to the thermal 
boosting case = 1.2, T+ = 1.1 for the Cased at = 1, k = 0.4, but are not sufficient 
enough to bring a significant change in load carrying capacity as seen in Tables 5.1 - 5.4. 



Figure 5.9: Pad temperatxires along the sliding direction of the bearing &t k = 0A,pf = 1 

In Figure 5.10 the u-velocity contours are plotted at = l,k — 0A,T^ = 1.3, — 1.3 

for Casel and Case2 respectively. Prom these Figures we see the manifestation of recircula- 
tion zones in region R. These zones can be attributed primarily to the thermal boundary 
conditions and partly to the geometric configurations of the slider bearings. 
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(a) (b) 


Figure 5.10: u-velocity contours at = 1, k = 0.4, = 1.3 for (a) Casel (b) Case2 

5.5 Conclusions 

Load carrying capacity and frictional drag of a high speed slider bearing with a solid pad 
for various thermal boundary conditions in conjunction with thermal boosting for various 
geometric configurations have been numerically analyzed using FEM. Density and the Vis- 
cosity in the mathematical model describing the slider bearing configuration are treated as 
functions of temperature. SUPG weight functions are found to be effective in precluding 
the numerical oscillations and in obtaining a convergent solution. Load carrying capacity of 
slider bearing can be thermally boosted by setting slider at temperatures lower than that 
of incoming lubricant. Though solid pad and slider conduction with different boundary set- 
tings as chosen in this study from application point of view alters the temperature field in 
the pad and thermal flux across the fluid-solid interface changes marginally and thus the 
load carrying capacity of the slider bearing remains unaffected. It is to be noted that the 
isothermal pad employed in the present study have been set at temperature similar to that 
of inlet flow. 



Chapter 6 


Thermohydrodynamic Analysis of a 
Tilted Pad Slider Bearing Considering 
Fluid Inertia 


6.1 Introduction 

In the previous Chapters, we have used Reynolds equation which explains the process of 
hydrodynamic lubrication in both qualitative and quantitative terms. The validity of the 
equation rests on the physical fact that the height of the fluid film is very small compared to 
its other physical dimensions. As a consequence all the hydrodynamic variations across the 
fluid film are negligible and the relevant Reynolds number in the bearings is usually small 
enough in case of laminar flows for the inertia and the turbulence effects to be neglected. 
At sufficient large Reynolds number and at high operating speeds and at low viscosity of 
the lubricant, the inertia force of the lubricant become of the same order of magnitude as 
the shearing force and the laminar flow gives way in stages to a turbulent state. However, 
even at moderate Reynolds numbers, for laminar flow of the lubricant, the contribution of 
the inertia force to the dynamics of bearings remain small. Thus, one visualizes (ref. Pinkus 
and Sternlich [74]) three broad regimes of the bearing operation: a lower regime of laminar 
flow in which the viscous forces are dominant (as we have seen in Chapters 3, 4 and 5), a 
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narrow intermediate regime in which the flow is still laminar but both viscous and inertia 
forces are significant, and an upper regime in which turbulent conditions prevail. The lower 
regime is exceptionally well described by Reynolds theory. On the other hand, in the narrow 
intermediate regime the existing theoretical formulations, like the iteration technique or the 
average inertia techniques, are applicable for one dimensional incompressible problems only 
and are not amenable to extensions. Apart from that these techniques cannot be looked upon 
as parts of a uniformly valid theory connecting the lower two regimes. As a consequence the 
utility of such scheme is very limited. 

In this Chapter, we have investigated the influence of convective inertia, various temper- 
ature boundary prescriptions and thermal boosting on bearing characteristics. The equations 
of mass, momentum and thermal energy conservations with temperature dependent density 
and viscosity .(which brings in extra pseudo compressibility terms into the energy equation) 
have been solved in a decoupled fashion using the SUPGFEM. 

6.2 Governing Equations 

=0 ( 6 - 1 ) 




du^ . + 0U+ , ^ 




^ ^ ^ 5^/+ ) P, dy+‘^ Pe dy+ p 

where, p+, h'^, k and P* are given by: 

p+ = 1 - A+ ( T+ - 1 ), = exp [ ( T+ - 1 ) ] 


dT+ , 1 d‘^T^ Pr Ea + . du+ .2 , Pr Ec + dp+ 


u 


dx+ 


(6.3) 


(6.4) 


= 1 — ( 1 — ^ ^ fi- — 


fio ^ T x>* pcJG" /k'o\2 




l^a B 


( 6 . 6 ) 
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The boundary conditions on the model are as follows: 

At inlet : u'*' = 1 - 

At outlet : 

On the pad : 

On the slider : 



= 0, 

p+ = 

Pi'^ 

(6.6) 





(6.7) 

dx'^ 

dy+ 


= 0 

u'^ 

■ = v^ = 

dp^ 

dfi 

= 0 

(6.8) 

= 

: I,!)"*" = 

dp+ 

dfi 

= 0 

(6.9) 


In view of the industrial applications, the energy equation has been solved under the fol- 
lowing different temperature boundary conditions: 

• Casel: (Isothermal slider and isothermal pad: ISIP or no flux) 

On slider: 0) =Ts^ 

On pad: T‘'’(x+, /i"*") = Tp^ 

In the current study, = Tp'*' = Ti^ 

• Case2: (Isothermal slider and adiabatic pad: ISAP or zero half flux) 

On slider: T+{x+, 0) = T/ (= Ti+) 

On pad: ^ (x+, /i+) = 0 

• CaseS: (Adiabatic slider and Adiabatic pad: ASAP or zero full flux) 

On slider: = 0 

On pad ^ = 0 

• Cased: (Isothermal slider and exposed pad: ISEP or non zero half flux) 

On slider: T+(x+, 0) = T/(= 7-+) 

On pad /c ^ = i? (T+ - 1)(= S) 

• CaseS: (Thermal boosting slider and pad: TBSP) 

(i). Tp+ = T/ > Ti+ (n).Tp+ > > Ti+ 

(iii).T/ > Tp+ > TZ {iv).Tp+=T+ < TZ 

The combination of T^, T+, which enhance the load carrying capacity of a slider bearing 
constitutes a thermal boosting model. At the inlet/outlet the following boundary conditions 
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are used: 


At inlet: T+(0, ?/■*■)= Tj 


+ 


At outlet: = 0 


The load carrying capacity (W'^) and the frictional drag force {F'^) have been determined 
by the following expressions: 

p+ dx'^ 

Jo 

du^ 

^ at 2 /+ = 0 


F+ = / ^ „+ _ A dx'^ 

J 0 


( 6 . 10 ) 

( 6 . 11 ) 


where, p"*", F'^ are given by: 

p* = l C = 


F+ = 


Fhp 
fJ^aU B 


6.3 Finite Element Formulation & Solution procedure 


Let fl C the domain of interest be bounded by a a piecewise smooth curve T. Let 
(a;"'’, y'^) denote any point in in the cartesian coordinate system and n be the outward unit 
normal drawn at any point (x'^,y'^) on F. Substitution for p'^, and into equations 


(6.2) and (6.3) and using equation (6.1) jdelds the following equations: 
[1-A+(T+-1)] (1^ + ^) - A+ (u+^ + 


'dx"^ 


[1 - A+(r+ - 1)] K C 


u 


dy+' 

_^. du+ 

dx"^ 


5x+ 


+ V 


+ 


dv'^ 

dy^ 


) = - 


dy+ 

dp^ 
dx+ 


(6.12) 


j. - B* 

^ ^ V ^ -1-2 ^ .4- .4- J 


[1 - A+(T+ - 1)] {u 


+ 


du^ 

dx'^ 


+ v' 


^dy+^ 
dv^ 


) = 


dy^ dy^ ‘ 
1 d'^T^ 


dy+^ Pe dy+2 


(6.13) 


+ 


Ec Pr 


nJPfm-P 


_L 


(6.14) 



6.3 Finite Element Formulation & Solution procedure 


119 


6.3.1 Petrov- Galerkin Formulation 


The Petrov-Galerkin weighted residual formulation of equations (6.12) - (6.14) gives the 
following equations: 






dp^ 

dx^ 




(6.16) 


£cPr j*(T+-i) ^ ^f2 ^ 0 

Pe da;+ 


(6.17) 


Pe ^ "dy^- 

Now to reduce the lubrication problem to one of finite number of unknowns, the solution 
domain is divided into bilinear rectangular elements and the field variables (pressure, veloci- 
ties and temperature) are expressed in terms of approximating trial functions Ni within each 
element. The domain now Q is discretized into Ne bilinear rectangular elements such that: 

u^i = (^ 


where, OS denotes the interior domain of an element. Let F® be the boundary of The 
discretized representation of the field variables are given by: 

Wei Net 

ut ^ Y. N.uf, vt « E W « E E’tpf, 

fc=l fc=l fc=l 

Wei 

Tt « E 

fc=i 


Considering the Petrov-Galerkin weighted residual formulation of equations (6.15) - (6.17) 
and introducing the above finite element discretizations we obtain the following: 
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/u.n. {II - - 1) ) IsIf ) 

^=1 k=l 

Nei Net N^I 

+ V - A+ 1 (E ■>4'm {-if (E I?" «■/=')) 

*=1 1=1 k=l 

N,i N,i 

+ (E (E } W'l rffj = 0 


(6.18) 


i=l 


k=:l 


iu.n- {(1 - A+(E (rr%) - 

k=l i=l fc=l 

ATe, Wei iVe, 

+ + idr&tPf) 


k=l 


k=l 


2=1 

iSTei Nel N^i 

expl-0*(ZTrNt - 

t=l k=l k=l 

Nel 

mda = o 

i=i 

Nel Nel Nel 

/u.n. {(1 - A+(E Wiv.) - i))((E«?‘JVf)5ff (Errjv« 

fc=l i=l ^=1 

Nel Nel Nel 

1=1 fc=i A=i 




- expl(-fi*(Y.^’N} - im^CZKO) 

i=l fc=l 

Nel Nel 

- w, df! = 0 

i=i fc=i 


where the Petrov-Galerkin weight function Wi is defined by: 

kujWij 

Iblf' 


^W^ = w^ + Pi,Wl = Ni,Pi 


(6.19) 


( 6 . 20 ) 


The upwind parameter k is calculated using elemental dimensions and elemental velocity 
(Grygiel and Tanguy [71], Brooks and Hughes [68]). Following the SUPG analysis of Brooks 
and Hughes [68], equations (6.19) and (6.20) can be cast into the weak form neglecting the 
inter element jumps arising out of definition of P,. This is justified as the bilinear elements 
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employed in the present analysis satisfy the required criteria like non-skewness as mentioned 
by Brooks and Hughes [68]. For solution evaluation a decoupled approach is taken. From 
equations (6.18) - (6.20) and boundary conditions the following decoupled equations are 
obtained for and T'^: 


(a) Decoupled equation for u"*" obtained from equation (6.16) is the following: 

^el fn-l) (n-1) 

E /u, n. { [(1 - A+(E ~ 1)) K ((E^.V ) 

k=l i=l 

^el (n-1) 

+ (E 

i=l 

^el ‘ (n— 1) (n — 1) 

+ 13* exp l-fi* (E rr Nf - 1)1 ^ (E N‘Tf ) ] W, 

1=1 j=l 

- exp - 1)] ^ p} dCl = 

i=l 

i=l 


/r w, exp 1-0* (E rr Ni‘ - 1)] ^ <ir 


(6.21) 


In view of the boundary conditions on the last term in equation (6.21) can be safely 
omitted without any loss of generality. In the above expression the superscripts (n-1) and 
(n) denote the values of the variables at the previous and current iterations and (n*) denotes 
the n level field variable values available at the time of the current calculation. 


(b) Equation for v'^ obtained from equation (6.18) is the following: 

Nel ^cl 


E fv.n. 


- JV,' ^ (E } W, vf' da 


- Su.A (1 - ^-"(E - 1)) sff (E ) 

j=l 

- AM (E m (p (E m ]} 

iz=l /c=l 


( 6 . 22 ) 
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(c) Equation for obtained from equation (6.20) is the following: 

^ , ,(«*) 


E /u.a. { (1 - ^-"(E Ni’Tf - 1)) ( (E 


k=l 


2=1 


2=1 


Net 




(n) 


+ (E ) 5? '^1 + T. W 


J = 1 




^el (n-1) pv- X 

/u,n. { ^ [-/?* (E - 1) 1 (E )) 

i=l 1=1 

t^el fn*) (n*) 

+ ^ (E (E )} 

i=l 1=1 

- Sr ^§WidT (6.23) 


Again in view of thermal boundary conditions excepting for the CaseS (ISEP) the boundary 
term on the right hand side of equation (6.21) may be omitted without any loss of generality, 
(d) Pressure is obtained from the momentum equation (6.19) which is analogous to Poisson 
equation for pressure from momentum equations. Equation for p'^ obtained from equation 


(6.19) is as follows: 


E /u.n. {Nt’Wi} pt <Sl 


k=l 


= - Su.u. { (1 - >^(E r^i’TF"~" -i))ir,w,{ (£N/uf “') 

i=l 1=1 

sff (E + (Eiv,'»r‘““’) ^iv.«(EAr,v'‘’“’) ) 

i=l 1=1 *=l 

+ expMME2^^"""V.^ - 1)] (^(E 0 

i=l 1=1 

+ P*w m ^ (EJV, %;•'”■”))} ■Ki 

i-l 1=1 

- /i,exp[-/3+(E - 1)1 (i)<^r 

2=1 


Nel 


(n-1) 


(6.24) 


Here again the boundary term on the right hand side of equation (6.24) may be safely omitted 
without any loss of generality. It must be mentioned that just to trigger the numerical 
calculations, pressures are initially obtained from the variational principle (Reddi [70]) for 
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Generalized Reynolds equation (Kumar et al. [73]). Now based on equations (6.21) - (6.24) 
the field variables are calculated in a decoupled fashion according to the 

following algorithm: 

6.3.2 Algorithm 

• Step: 0, Initialization: 

(a) Carry out boundary settings for 

(b) Set initial lubricant temperature to inlet lubricant temperature. 


Step: 1, If i = 1 (first iteration) then Evaluate p'*’ minimizing the variational principle 
for Generalized Reynolds equation (Kumar et al. [73]). i.e. 

<5n(p) = / [-^ ( 1 ^) - ^] dQ (6.25) 

^ Ja 12p+ ^dx+ dx+^ 


else 

Evaluate based on equation (6.24) using 




Step: 2, If i = 1 (first iteration) then 

Evaluate using from the direct integration of linearized x-momentum equa- 
tion i.e. 


d , + dw^ . _ dp'*’ 
dy+ ^ ~ dx+ 


(6.26) 


else 

Evaluate evaluate using 


• Step: 3, If i = 1 (first iteration) then 

Evaluate directly from the continuity equation treating p'*’ to be a constant and 
using 

else 

Evaluate based on equation (6.22) using \ 
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The convex parabolic pressure profiles corresponding to {pt,pt) = (0, 0) clearly indicate 
the pressure generation in the slider bearings for zero inlet pressure. It can be easily observed 
that the increase in the tilt of the pad significantly enhance the pad pressures and shifts the 
location of peak pad pressure occurrence closer to the outlet. These observations are very 
much in tune with the earlier studies (Pinkus and Strenlicht [74]). To study the effect of 
non zero inlet pressure, simulations are carried out maintaining (pt,Po) at (1,0) and fixing 
2t+ = T^ — = 1.3 for /c = 0.4, 0.6, 0.8, 1.0. Pad pressures along the sliding direction are 

shown in Figure 6.2. 



Figure 6.2: Pressures along the pad for pf = 1,T^ = = 1-3 

Again the pad pressure increase with the tilt in the pad and the peak pad pressures are 
noticed closer to the outlet. For k = 0.4, 0.6, 0.8 it is to be noted that the peak pressures 
seen near the outlet are higher than the prescribed inlet pressure. These observations are 
compatible with several experimental findings Pinkus [31] and indicate the possibilities of 
enhancing load carrying capacity of slider bearings by considering fore-region pressure and 
necessitates further investigations on the influence of fore-region pressure. Hence simulations 
are carried out for pf = 2.5,5,10 fixing Tf' = Tf' = Tf' = 1.3 for k = 0.4, 0.6, 0.8, 1.0. 
In Figure 6.3 pressure distribution on the pad in sliding direction for k = 0.4 and pf = 
0, 1, 2.5, 5, 10 are presented. 
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Figure 6.3: Pad pressures k = 0.4, =T^ = T+ = 1.3 for various values of 


Clearly an increase in the inlet lubricant pressure increases the pad pressures which in 
turn leads to an enhanced load carrying capacity. The load carrying capacity of the slider 
bearing corresponding to the above settings are evaluated and are presented in Figure 6.4. 



Figure 6.4: Influence of various inlet pressures on load carrying capacity for = T+ = 1.3 

for various values of k 

From this Figure 6.4 it may be noted that the load carrying capacity of slider bearing 
can be enhanced either by increasing the tilt of the pad (fc) or by the consideration of fore- 
region pressure {pf ) . Clearly the combined influence of k and pf on the load generation is 
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greater than their segregated influences. 

The frictional drag coefiicient associated with above setting are evaluated and presented 
in Figure 6.5. 



Figure 6.5: Influence of various inlet pressures on frictional drag force for id" = y + = T+ = 1.3 
for various values of k 

It is amply clear that the frictional drag reduces with the tilt in pad and increases 
with the increasing This may be attributed to the acceleration in the flow brought in 
by fore-region pressures or by the geometric effect of the slider bearing. All of the above 
investigations pertain to the ISIP settings. It is essential to investigate the influence of other 
boundary settings as defined under cases ISAP, ISEP and ASAP on load generation and 
frictional drag. The load carrying capacity and the frictional drag associated with the slider 
bearing for these various boundary settings are evaluated and presented in Tables 6.1 and 
6 . 2 . 

Here the computations have been done for {pf ,p^) = (0, 0), (1, 0). ,T ^ ,Tj are fixed 
at 1.3, whenever applicable and on ad-hoc out flux of magnitude 0.2 is considered for ISEP 
(non zero half flux) case. From Table 6.1 one can notice that the consideration of fore-region 
pressures brings in a significant enhancement in load carrying capacity for all boundary 
settings. Further the boundary settings corresponding to the case ISEP (non zero half flux) 
lead to relatively higher load carrying capacity of the slider bearing. 
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Table 6.1: Load capacity at various settings for = T/ = 1.3 


Pi 

w+ 

w+ 

w+ 

w+ 

Cases 

ISIP 

ISAP 

ASAP 

ISEP 

k = 1 

0.0003775 

0.00037653 

0.0037654 

0.000399943 

k = 0.8 

0.007463 

0.007456 

0.007455 

0.0094824 

II 

p 

0.0188987 

0.01889406 

0.0188956 

0.019489 

k = 0.4 

0.04600167 

0.0460329 

0.0460306816 

0.04702043 


w+ 

w+ 

w+ 

w+ 

k = 1 

0.559023 

0.5590082 

0.5590084 

0.5751408 

k = 0.8 

0.69404 

0.69621 

0.69622 

0.698157 

k = 0.6 

0.741607726 

0.74160432 

0.741603315 

0.7423247 

o 

II 

0.8276344 

0.82726675 

0.8272555 

0.8539077 


Table 6.2: Frictional drag force at various settings for T+ = = T/ = 1.3 


pt = 0 

F+ 

F+ 

F+ 

F+ 

Cases 

ISIP 

ISAP 

ASAP 

ISEP 

k = 1 

-0.049685374 

-0.049685359 

-0.0496150441 

-0.049724504 

k = 0.8 

-0.058889322 

-0.058889363 

-0.0587618388 

-0.0588530414 

k = 0.6 

-0.0746395215 

-0.074275121 

-0.074275121 

-0.0746286213 

k = 0.4 

-0.103227668 

-0.103228159 

-0.102606408 

-0.103196517 

pt = 1 

F+ 

F+ 

F+ 

F+ 

k= 1 

-0.03409398 

-0.0340939611 

-0.034287925 

-0.0341370404 

k = 0.8 

-0.0394585505 

-0.039458595 

-0.038821835 

-0.0392912589 

k = 0.6 

-0.0609837808 


-0.0614150167 

-0.060973309 

k = 0.4 

-0.0913590118 

-0.0913589522 

-0.0907316506 

-0.0913887247 


The load carrying capacity of slider bearing do not change or differ marginally for 
ISEP, ISAP and ASAP boundary settings. Recall that the viscosity and the density of the 
lubricant given by equation (6.4) are sensitive to the temperature. From the expression for 
the viscosity it is clear that for fixed values of /3’*' smaller the value of (T”^ — 1) larger vpill 
be Generally, lubricants with larger viscosity and density will have higher threshold for 
thermal expansion and can lead to relatively larger load carrying capacity. So the lubricant 
temperature fields associated with the ISEP, ISIP, ISAP and ASAP boundary settings can 
explain the observed behavior of the load carrying capacity of slider bearings as one moves 
from ISIP (Case 1) settings to ISEP (Case 4) settings. Hence in Figures 6.6 and 6.7 lubricant 
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From the temperature fields as seen in Figures 6.6 and 6.7 one can notice that the 
average temperature of the lubricant for ISIP, ISAP, ASAP cases is above 1.3. Whereas for 
ISEP it is around 1.2. Further from the Figures 6.6 and 6.7 corresponding to ISIP, ISAP 
and ASAP cases, hot temperature zones where the temperatures rises above 1.3 are noticed 
in the regions 0.2 < a; < 0.4, 0.3 < y < 0.6 and 0 < a: < 0.5, 0.7 < y < 1.0. But in 
Figures 6.6 (d) and 6.7 (d) corresponding to ISEP boundary settings no such hot spots are 
noticed. Higher temperatures lead to lower effective viscosity and density. Consequently 
one gets a lower load carrying capacity of the slider bearing associated with the ISEP, ISIP, 
ISAP boundary settings. The temperature distribution in the Figures 6.6 and 6.7 suggests 
that ISEP boundary settings can lead to better performance of slider bearing. To ensure 
that these conclusions and reasoning independent of k, isotherms for is = 1 at pf = 0, 1 
for the boundary settings ISIP, ISAP, ASAP and ISEP are presented in Figures 6.8 and 6.9 
respectively. Again observations similar to those found in Figures 6.6, 6.7 are noticed. The 
circular isotherms found in the above plots indicate the possibility of recirculation zones in 
the interior region of flow field. 

To further investigate these recirculation zones the iso-u velocities or streamlines (as 
v-velocities are negligible) for k = 0.8, 0.4 at pf = 0, 1 respectively are presented in Figures 
6.10 and 6.11. From these Figures it may be seen that for all the four boundary settings and 
both at = 0, 1 circular contours indicating recirculation zones are noticed in the region 
0.2 < a; < 0.4, 0.3 < y < 0.6. Exactly this is the region of high temperature zone or hot 
spot in the corresponding lubricant temperature field. The combined influence of geometric 
parameter k and the consideration of fore-region pressure on the flow field is clearly visible 
as the laminar boundary aligned streamline flow seen at = 0 for k — 0.8 changes to a 
gushing flow at pf^ = 1 for k = 0.4. This also suggests that for high speed lubrication one 
can not rely on conventional models based on linearized velocity calculations. 






Figure 6.9: Iso-temperature lines for pf — I, k — 1-0, 
(c) Cases and (d) Case4 


= 7 + = T+ = 1.3; (a) Casel (b) Case2 
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The frictional drag coefficient associated with ISIP, ISAP, ASAP and ISEP boundary 
settings for k = 0.4, 0.6, 0.8, 1.0 at pf = 0, 1 (as seen in Table 6.2) depict only a marginal 
variation. Recall that the frictional drag co-efficient is evaluated using equation (6.11) and 
viscosity is temperature dependent. Now a look at the flow field and temperature field in 
the immediate vicinity of y'^ = 0, where frictional drag co-efficient is calculated, show hardly 
any variation, as one moves from ISIP setting to ISEP settings and hence the marginal 
differences in the frictional co-efficient. 

For high speed slider bearings inclusion of inertia effects, as in the present case, has 
been recommended by several researchers (Rodkiewicz and Sinha [58]). To understand the 
influence of inertia effects on load carrying capacity a comparison of load carrying capacity 
obtained in the current simulation wherein inertial effects have been considered with the 
earlier simulations devoid of inertial effects is made in Table 6.3. 

Table 6.3: Comparison of load capacity with and without fluid inertia 


p+ = 0 





k 

1 

0.8 

0.6 

0.4 

W'^ (without inertia) 

0.000000464727 

0.004307854 

0.012646012 

0.03063336 

W'^ (with inertia) 

0.00003774193 

0.006632 

0.0188987 

0.04600167 

pt = i 


W'^ (without inertia) 

0.4999 

0.5580224 

0.63334724 

0.7378243 

W'^ (with inertia) 

0.559023 

0.5702698 

0.6403325 

0.8539077 


It is interesting to note that consideration of inertial effects enhance the load carrying 
capacity of high speed slider bearings for all values of k and 

The influence of thermo viscosity parameter {j3'^) on the load generation in high speed 
slider bearing lubrication has been analyzed at different fore-region pressures for various 
slider bearing geometries. Plots in Figure 6.12 depict the dependency of load carrying 
capacity on /J'*' for = 0,0.5, 1. For small fore-region pressure load generation in high 
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speed slider bearings decrease with increasing /?+ but at high inlet pressures this reduction 
in load generation becomes marginal. This may be explained directly from the expression 
for ij,'^ given in equation (6.4). It can be noted that for fixed T+,y^+(T+) increases with 
decreasing 0'^, and thus the effective viscosity of the lubricant increases and will lead to 
higher load carrying capacity. 



Figure 6.12: Comparison of load carrying capacity with yS"*" at = 1.3 for various 

values of k and pf 

Now to find the thermal boosting boundary settings for slider bearing simulations have 
been carried out at pf = 1, fixing (a) = T+ = 1.3, = 1.2 (b) Jd” = T/ = 1.3, T+ = 

1.2 (c) T+ = T+ = 1 . 2 , 27 ^ = 1.3 (d) = T+ = T+ = 1.3 and are shown in Table 6.4. 

Table 6.4: Load capacity at various settings for pf = l,k = 0.4 


"jiT' 

-*•5 

rp-j- 

w+ 

1.3 

1.2 

1.3 

0.795805991 

1.3 

1.3 

1.3 

0.85390777 

1.3 

1.3 

1.2 

1 S 

1.3 

1.2 

1.2 

0.887764513 


From the Table 6.4 we see that the thermal settings under (b) and (c) can boost load 
carrying capacity of high speed slider bearings. Maintaining either pad at lower temperature 
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or both pad and slider at lower temperatures than the incoming lubricant thermally boost 
the load carrying capacity of the slider bearing. To obtain a causative picture of this thermal 
boosting of load generation, temperature fields in the form of isotherms associated with the 
above settings are presented in Figures 6.13 and 6.14. From the temperature distribution 
seen in these Figures, the net average lubricant temperatures associated with the above cases 
(a), (b), (c) are in the order Tf^a) > T(b) > T(^c) > 1- Now by the reasoning presented earlier 
i.e. based on the trend observed in the load carrying capacity gets justified. It may 

be noted that the pad is stationary and it is reasonable to find that maintaining it at lower 
temperature leads to increased load generation. 



(a) (b) 


Figure 6.13: Iso-temperature lines for = 1.3,/: = 0.4,p^ = 1, Casel: (a) T+ — 1.3, — 1.2 

(b) r+ = i.2,r+ = i.3 
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Figure 6.14: Iso-temperature lines for = 1.3, A; = 0.4, = 1, Casel: T+ = T+ = 1.2 


6.5 Conclusions 

The thermal effects and inertia effects on bearing characteristics in slider bearings have been 
analyzed. An orderly analysis identifying the influence of various factors viz. inlet pressure, 
inlet lubricant temperature, thermal boosting, slider bearing geometry on load carrying 
capacity and frictional drag force have been provided. The Figures and Tables depicted 
in this Chapter reveal that the consideration of inlet pressure or thermal boosting or both 
will enhance the load carrying capacity of a shder bearing with a reduction in frictional 
drag force. Further, it has been observed that consideration of fluid inertia will enhance the 
effective load carrying capacity of a slider bearing. 



Chapter 7 


Thermohydro dynamic Analysis of a 
Tilted Pad Slider Bearing Considering 
Fluid Inertia and Heat Conduction in 
the Pad 


7. 1 Introduction 

It has already been mentioned in Chapter 4 and also in Chapter 6 , that for lubrication 
under high loads and high speeds, large amount of heat may be generated. Heat created in 
the lubricant film has a profound effect on bearing characteristics. Thus, an analysis, aimed 
at accurately predicting the bearing performance must consider not only the conservation 
equations of mass, momentum in the fluid film, but also the conservation equation of thermal 
energy in the fluid as well as conduction in the stationary pad. 

Studying the nature of thermal effects and heat conduction in the stationary pad, Stren- 
licht [75] performed a restricted adiabatic analysis of a sector-shaped thrust bearing. His 
results showed that temperature variations in the direction perpendicular to the direction of 
sliding are not always negligible. 

Although studies up to 1957 did not allow for temperature variations across the film, 
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they did demonstrate that accurate prediction of bearing performance is seldom possible 
if fluid-film property variations are neglected. Zienkiewicz [42] was the first to formulate a 
numerical scheme for simultaneously solving the conservation equations in the fluid film with 
temperature variations across the film. Zienkiewicz’s and later follow-up results of Hunter 
and Zienkiewicz [12] convincingly indicate that viscosity variation with temperature is far 
more important than density variation. 

In a pair of papers, Dowson and Hudson [44] performed a thermohydrodynamic analysis 
of parallel and inclined slider bearings of infinite length. Their results dramatically show 
the negligible effect of density changes and the dominant effect of viscosity changes due to 
temperature gradients across the film. 

Hahn and Kettleborough [45] analytically approached the “Fogg Paradox” (ref. Fogg [36] 
by considering thermal distortion and viscosity variations along and across the film. They 
have presented an analysis for heat conduction and thermal gradients in one-dimensional 
slider bearings having a general film shape. 

In this Chapter, we have analyzed the slider bearing performance by considering the 
governing equations of mass, momentum and energy in the lubricant film and the conduction 
equation in the solid pad, using the SUPG Finite Element Method. 

7.2 Governing Equations and Solution Procedure 


K p* ( « 


+ 


dw^ 

dx+ 


+ v"' 


du'^ 

dy+ 


) = 


dp 


,+ 




+ 


d 


,+du+ 


dy+ dy+ 


(7.2) 




dT+ 

dx^ 


-f v'^ 


dT+ _ J_ 
dy+ ^ ~ Pe 


PrEc + . du+ 
^ Pe “ ^ 




dp'^ 

da:+ 


(7.3) 
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where, A: and are given by: 




1 - A+ ( T+ - 1 ), /i+ = exp [ ( r+ - 1 ) 


(7.4) 


h* = l-x* (l-k), k=^ < 1, Rl = 

hi B 


(7.5) 


The equation governing heat conduction in the stationary pad is: 

rP'T + SP'T + 

+ =0 
dx'^ dy^ 


(7.6) 


Now using the Petrov-Galerkin finite element discretization for the above equations, we 


obtain the following equations. 


ST+ 5T+ 


/ {[i-A+(T+-i)jij:(«+ 

«/ 0 




Sr;'*’, 


dp 


,+ 


^ Sy*" da:’*’ 


-^+(r+-i) 


^Sy*’^ 


- 




Sy+ Sr/’*’ 


)}T^dfl = 0 


Sr;’*'. 


J_ S^T+ 
Sy*" Pe Sy*'^ 


(7.8) 




_ ScPr ^ <jn = 0 


Pe 


Pe 


da:+ 


(7.9) 


And the Galerkin finite element discretization for heat conduction equation is as follows: 

L ' ^ ^ ay 

Equations (7.7)- (7. 10) are solved iteratively (for details see Chapter 6) under the following 

temperature boundary conditions (see Figure 3.1), to yield the field variables viz. ■u'*', r;’*’, p'*’, T’*'. 
• Casel: 


Lubricant boundary conditions: 
T+ = T+ = 1.3, on AB 
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dT+ 


0, on BC 


T+ = T^ = 1.3, on AD 
Pad boundary conditions: 

Tp+ = 1.3„ onDF 
T+ = 1.3„ onCE 
T+ = 1.3„ onEF 
Interface boundary conditions: 

T"*" = Tp'^ and kg = kf on 
• Case2: 

Lubricant boundary conditions: 

T+ = T+ = 1.3, on AB 


^=0, onBC 
r+ = I^+ = 1.3, on AD 
Pad boundary conditions: 

^ = 0, onDF 
^ = 0, CE 

dx'^ ’ 

Tp+ = 1.3, onEF 
Interface boundary conditions: 

T+ = T* and K^ = k 
• Case3: 

Lubricant boundary conditions: 
T+ = T+ = 1.3, on AB 


/ dy^ ’ 


on 


fj = 0. onBC 

=Tt = 1.3, on AD 
Pad boundary conditions: 
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^ = 0, onEF 
Interface boundary conditions: 

T+=r/ and on CD 

• Cased: 

Lubricant boundary conditions: 

T+ = T+ = 1.3, onAB 
iJ = 0. onBC 
T+ = = 1.3, on AD 

Pad boundary conditions: 

1^ = 0, onDF 
^ = 0, onCE 

= (V-1), onEF 

Interface boundary conditions: 
r+ = T/ and onCD 

7.3 Results and Discussion 

The simulations have been carried out on 30 x 30 grid system. The results have been 
presented for load carrying capacity in the form of tables for various boundary settings. 

The load carrying capacity for the setting = T/ = 1.3,p^ = 1 have been 

compared with the results obtained for a linearized velocities model (Chapter 4) and are 
shown in Table 7.1. Inclusion of inertia results in an enhancement of load capacity for all 
values of k. 

In Table 7.2, comparison of load carrying capacity with conduction in the pad (Chapter 
7) and with no heat conduction in the pad (Chapter 6) have been made for various slider 
bearing geometries with the setting = T/ = 1.3,p^ = 1 Case4 boundary condition. 
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Table 7.1: Comparison of load capacity with inertia and without inertia forces for pf = 1, Case4 


Cases 

II 

O 

k = 0.6 

k = 0.8 

k = 1.0 

W"*", without inertia (linearized model) 
W'^, with inertia (nonlinear model) 

0.7378243 

0.855990 

0.633222824 

0.7502945 

0.55832 

0.69157 

0.49989 

0.57481 


Table 7.2: Comparison of load capacity with conduction and without conduction in the pad for 
= 1, Case4 


Cases 

k = 0.4 

CD 

II 

k = 0.8 

k = 1.0 

with conduction 
W*" without conduction 

0.855990 

0.8539077 

0.7502945 

0.7423247 

0.691857 

0.698157 

0.57481 

0.575108 


Prom Table 7.2, It is observed that the load carrying capacity associated without con- 
duction in the pad are marginally higher for the cases k = 1,0.8, whereas for the case 
k — 0.6, 0.4, these are lower. It shows that as the tilt of the pad increases the load carrying 
capacity decreases for the case with conduction in compared with the case without conduc- 
tion. Similar observations were made for the linearized model ( Chapters 3, 4). Inclusion of 
the inertia forces and heat conduction in the pad supports the pressure generation when the 
tilt is raised from k = 0.6 to 0.8. 

In Table 7.3, the load carrying capacity for various boundary conditions at various values 
of k for the setting = T+ = 1.3, = 1 are shown. From Table 7.3, one can observe 

that Cased boundary condition has higher values of load carrying capacity for various values 
of k as compared with the other Cases. This is due to the fact that the temperature of 
the lubricant decreases when we allow the pad conduction to the environment. Also, it is 
observed that for parallel slider case (k = 1) there are negligible variations in the computed 
load carrying capacity for all the temperature boundary conditions. This could be due to 
no change in the temperature gradients as we have seen in Chapter 5. 
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Table 7.3: Comparison of load capacity for various boundary conditions pf = = 

T+ = 1.3 


Cases 

k = 0.4 

oo 

o 

II 

II 

p 

O 

t-H 

II 

W+, Casel 

0.8350990 

0.735544 

0.6947003 

0.559120 

W+, Case2 

0.824330 

0.725974 

0.6947007 

0.559001 

W+, Cases 

0.8251006 

0.725946 

0.690537 

0.559008 

W+, Case4 

0.855990 

0.7502945 

0.691857 

0.575481 


7.4 Conclusions 

Load carrying capacity of a high speed slider bearing with conduction in the stationary 
pad for various thermal boundary conditions for various geometric configurations have been 
numerically analyzed using FEM. Density and the Viscosity in the mathematical model 
describing the slider bearing configuration are treated as functions of temperature. Though 
solid pads with different boundary settings as chosen in this study from application point of 
view alters the temperature field in the pad the thermal flux across the fluid-solid interface 
changes marginally and thus there is a little change in the load carrying capacity of the slider 
bearing. 
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